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Abstract 

The  objective  of  the  current  research  was  to  investigate 
the  effects  of  reducing  screen  thickness  on  the  volume  and 
compactness  factor  ijn/f)  of  stacked,  wire-screen 
regenerators.  An  improved  transient  step-change  method  was 
devised  which  integrated  experimental  data  with  a  numerical 
model  of  the  flow  to  determine  the  heat  transfer  coefficient 
and  friction  factor.  The  improvements  to  the  approach  are: 

1)  the  measured  inlet  temperature  trace  is  used,  2)  the  heat 
transfer  coefficient  is  based  on  a  parameter  called  the  sponge 
effect  delay  time,  and  3)  the  important  effect  of  the  tube 
surrounding  the  matrix  is  included  in  the  numerical  model. 
The  data  show  that  the  heat  transfer  is  the  same  for  reduced 
thickness  screens  as  it  is  for  unrolled  screens  once  the 
decrease  in  surface  area  caused  by  rolling  the  screens  is 
taken  into  account.  However,  the  friction  factor  increases, 
particularly  for  a  50%  reduction  in  screen  thickness. 
Consequently,  the  ratio  of  Colburn  factor  to  friction  factor, 
the  compactness  factor,  decreases  as  the  thickness  of  the 
screens  decrease.  The  effectiveness  of  the  regenerators  is 
also  adversely  affected  by  the  rolling. 


AN  INVESTIGATION  OF  THE  CHARACTERISTICS 


OF  REGENERATIVE  HEAT  EXCHANGERS 

X-  Introduction 

As  the  chapter  title  states,  this  is  an  introduction  to 
the  research  described  in  this  dissertation.  A  concise 
statement  of  the  problem,  motivation,  and  approach  for  the 
research  is  given  here,  as  well  as  two  definitions  which 
require  clarification  in  advance.  Background  information  on 
stacked-screen  regenerators  and  justification  for  the  research 
are  given  in  Chapter  II.  Chapter  III  describes  the  theory 
behind  determining  a  heat  transfer  coefficient  for  a  porous 
medium,  and  describes  the  improved  approach  used  for  the 
current  research.  The  experimental  apparatus  and  procedure 
are  described  in  Chapter  IV.  The  results  are  presented  in 
Chapter  V,  while  some  concluding  remarks  and  recommendations 
appear  in  Chapter  VI .  Other  important  information  which 
includes  MATLAB  m. files,  FORTRAN  codes,  test  procedures, 
calibration  techniques,  and  experimental  data  are  given  in  the 
appendices . 

Obi ective 

The  hypothesis  central  to  this  research  may  be  stated  as 


follows : 


The  reduced  thickness  of  pressed,  round  wire 
screens  decreases  the  volume  of  a  stacked, 
wire-screen  regenerator,  and  reduces  its 
compactness  factor  (Jh//) • 

The  definition  of  pressed  is  any  process  which  reduces  the 
thickness  of  the  screens  without  removing  material.  The 
specific  objective  is  to  determine  the  heat  transfer 
coefficient,  friction  factor,  and  other  flow  and  geometrical 
properties  of  a  porous  medium  regenerator  made  from  rolled 
screens,  and  to  compare  these  results  to  those  for  unrolled 
screens.  To  accomplish  this  objective,  an  integrated 
experimental-numerical  method  which  produces  dependable  data 
was  deve 1 oped . 

Motivation 

Improving  the  performance  of  regenerative  refrigeration 
cycles  is  important  to  the  Air  Force.  A  better  understanding 
of  the  processes  in  porous  medium  regenerators  which  are  an 
integral  part  of  regenerative  refrigeration  cycles  will  lead 
to  improvements  in  performance.  One  specific  improvement  is 
to  reduce  the  dead  volume  of  the  regenerator.  But  this  should 
only  be  done  if  the  method  for  reducing  the  dead  volume  does 
not  decrease  the  heat  transfer  coefficient  or  increase  the 
pressure  drop  characteristics  of  the  regenerator 
significantly,  since  good  performance  also  depends  on  good 
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heat  transfer  and  small  pressure  drops.  Hence,  a 
comprehensive  study  of  the  effects  of  reducing  the  screen 
thickness  on  the  compactness  factor,  which  is  the  ratio  of  the 
Colburn  factor  to  the  friction  factor  (Radebaugh  and  Louie, 
1986:180),  was  conducted. 

Approach 

A  three-step  approach  to  this  research  was  used.  First, 
an  experimental  apparatus  was  designed  and  built  to  measure 
temperatures,  pressures,  and  other  properties  of  reduced¬ 
thickness,  wire-screen  regenerators  for  a  range  of  flow  and 
geometrical  conditions.  Second,  a  data  reduction  technique 
and  numerical  model  of  the  flow  that  determine  the  heat 
transfer  coefficient  based  on  criteria  from  the  experimental 
data,  and  that  calculate  friction  factors  for  each  case,  was 
developed.  And  third,  the  results  were  analyzed  and 
conclusions  about  the  effect  on  the  compactness  factor  of 
reducing  the  volume  by  rolling  the  screens  were  drawn.  A  by¬ 
product  of  this  approach  was  a  better  understanding  of  how  the 
heat  transfer  coefficient  can  be  more  accurately  determined  by 
accounting  for  non-idealizations  and  omissions  present  in  the 
experimental  and  analytical  techniques  used  by  other 
researchers  of  regenerative  heat  exchanger  characteristics . 
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Important  Terms 


Two  terms  should  be  clarified  in  advance.  First,  the 
matrix  is  the  porous  medium,  i.e.  the  stack  of  wire-mesh 
screens.  The  regenerator  is  a  collection  of  parts  which 
includes  the  matrix,  tube,  spacers,  connectors,  etc.  Second, 
SU  is  the  designation  for  screen  unit.  Thirteen  screen  units, 
SUs,  were  fabricated  for  the  current  research,  each  with  a 
different  combination  of  mesh  size  and  screen  thickness 
reduction  factor.  More  definitions  of  terms  are  made  in 
Chapter  II,  Background  and  Justification. 
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II .  Background  and  Justification 


The  current  research  is  being  done  to  contribute  to  the 
knowledge  of  heat  and  momentum  transfer  in  a  porous  medium, 
and  to  demonstrate  a  better  way  to  study  these  phenomena.  The 
how  part  will  be  addressed  in  Chapter  III,  Methodology.  This 
chapter  describes  the  rationale  for  the  hypothesis  stated  in 
Chapter  I .  The  Background  section  relates  the  importance  of 
cryocooler  technology  to  the  Air  Force,  defines  terms,  and 
explains  how  regenerators  are  essential  to  the  operation  of 
cryocoolers .  In  the  Justification  section,  the  literature 
search  and  available  data  are  described,  and  the  hypothesis  is 
formulated. 

Background 

The  particular  application  which  motivates  the  current 
research  is  regenerative  refrigeration.  Regenerative 
refrigeration  technology  is  important  for  a  wide  range  of 
space-based  applications  (Chan  et  al,  1990;  Thieme  and  Swec, 
1992;  Ledford,  1994).  The  Air  Force  currently  oversees  a 
large  research  and  development  program  in  this  area  (Thomas, 
1992;  Ross,  1992;  Wyche  and  Bruning,  1990).  Regenerative 
refrigeration  cycles  are  especially  suited  to  space 
applications  because  they  are  low-maintenance,  low-vibration 
devices  with  demonstrated  long-endurance  capabilities  (Walker, 
1983:  Chapter  4) .  At  the  low  temperatures  commonly  found  in 
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outer  space  (less  than  120  K)  ,  these  refrigerators  are 
referred  to  as  cryocoolers  (Radebaugh,  1989)  .  The  most 
important  application  for  cryocoolers  is  in  thermal  management 
of  electronic  components.  In  particular,  satellite  imaging 
systems  are  equipped  with  semi-conductor,  photon-collecting 
components  that  are  kept  as  cold  as  possible  for  maximum 
resolution  (Chan  et  al,  1990) . 

The  stringent  requirements  placed  on  space-based,  energy 
conversion  equipment  exist  for  one  simple  reason:  Routine 
maintenance  and  refueling  are  impractical  for  space-based 
assets  due  to  a  lack  of  access.  Consequently,  even  small 
improvements  in  energy  conversion  efficiency  are  important  for 
extending  space  asset  lifetimes.  Much  of  the  current  research 
on  cryocoolers  pursues  this  goal  (Thieme  and  Swec,  1992). 

There  are  various  types  of  regenerative  refrigeration 
cycles  used  in  cryocoolers,  but  the  free-piston  Stirling  cycle 
has  received  the  most  attention  since  it  is  very  reliable  and 
holds  promise  for  high-efficiency  space  utilization  (Chan  et 
al,  1990:  1244;  Ledford,  1994).  As  the  name  regenerative 
refrigeration  cycle  implies,  a  regenerator  is  an  integral  part 
of  the  system.  In  order  to  demonstrate  the  purpose  and 
importance  of  the  regenerator,  a  Stirling  refrigeration  cycle 
will  be  examined. 

The  Stirling  Cycle.  This  section  describes  the  basic 
elements  of  a  Stirling  refrigerator  and  the  associated  ideal 
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thermodynamic  cycle.  The  aim  of  this  section  is  to  explain 
what  a  cryocooler  is  and  how  it  goes  about  performing  its 
task . 

Fig.  1  shows  the  basic  elements  of  a  simple,  two-cylinder 
Stirling  cycle  refrigerator.  The  working  fluid  is  typically 
a  low  boiling-point,  monatomic  gas  such  as  helium.  Heat  is 
transferred  from  the  refrigerated  volume  at  Tg  into  the 
expansion  space;  heat  is  rejected  to  a  sink  at  T^  from  the 
compression  space.  A  regenerator  is  placed  between  these  two 
spaces.  The  gas  is  under  high  pressure,  typically  around  25 
atmospheres,  and  oscillates  between  the  compression  space  and 
the  expansion  space  through  the  regenerator.  A  regenerator  is 
a  heat  exchanger  which  receives  energy  from  the  gas  during  one 
part  of  the  cycle  (the  hot  blow)  ,  holds  the  energy,  and 
returns  it  to  the  gas  on  the  reverse  passage  through  the 
regenerator  (the  cold  blow) .  The  regenerator  may  be  thought 
of  as  a  thermal  sponge,  alternately  absorbing  and  releasing 
heat  energy  (Walker,  1983:45). 

The  thermodynamic  processes  of  an  ideal  Stirling  cycle 
are  shown  in  the  P-V  and  T-S  diagrams  in  Fig.  2.  The  cycle 
begins  with  all  of  the  gas  in  the  compression  space.  A 
constant-temperature  compression  occurs  between  state  1  and 
state  2  during  which  heat  is  transferred  out  of  the  gas  at  T,,. 
Next,  a  constant-volume  heat  transfer  takes  place  during  which 
heat  is  transferred  from  the  gas  to  the  regenerator.  From 
state  3  to  state  4,  a  constant-temperature  expansion  occurs 
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during  which  heat  is  transferred  into  the  gas  from  the 
refrigerated  volume  at  T^.  Finally,  the  energy  stored  in  the 
regenerator  during  the  hot  blow  is  returned  to  the  gas  by  a 
constant -volume  process  which  returns  it  to  state  1. 
Differences  exist  between  real  and  ideal  Stirling  cycles,  but 
the  above  description  is  sufficient  for  illustrative  purposes. 

The  dominant  energy  transfer  process  in  a  Stirling  cycle 
is  the  exchange  of  heat  between  the  working  gas  and  the 
regenerator.  The  heat  transferred  to  the  regenerator  during 
each. half -cycle  is  typically  ten  to  fifty  times  larger  than 
the  heat  removed  from  the  refrigerated  volume  (Atrey  et  al .  , 
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Figure  2.  Ideal  Stirling  Thermodynamic  Cycle 

(Walker,  1983:45) 


1990:236).  The  heat  transferred  to  the  regenerator  during 
each  half-cycle  is  usually  even  larger  than  the  work  input 
from  the  compressor.  Consequently,  the  ability  of  the 
regenerator  to  effectively  transfer  and  store  energy  is  veiry 
important  to  the  performance  of  a  cryocooler. 

Cryocooler  performance  can  be  measured  in  various  ways . 
A  common  figure  of  merit  for  a  refrigeration  cycle  is  the 
coefficient  of  performance  (COP) .  For  a  refrigeration  cycle, 
work  is  put  into  the  compression  space  and  heat  is  removed 
from  the  refrigerated  volume.  The  COP  for  a  cycle  is  defined 
as  the  ratio  of  the  heat  removed  from  the  refrigerated  volume 
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to  the  net  work  input  to  the  compression  space  (Walker, 
1983:42) .  Consequently,  a  high  COP  is  desirable. 

The  COP  is  a  system-level  f  igure-of-merit  for  a 
cryocooler.  For  the  regenerator  component  of  the  cryocooler, 
the  compactness  factor  is  used  as  the  f igure-of-merit .  It  is 
a  ratio  of  the  heat  transfer  to  pressure  drop  in  the 
regenerator.  The  heat  transfer  is  characterized  by  the 
Colburn  factor,  jjj,  while  the  pressure  drop  is  characterized 
by  the  friction  factor,  f  (Kays  and  London,  1984:7). 

A  relationship  between  the  COP  and  the  compactness  factor 
can  be  written  for  a  given  set  of  operating  conditions  and 
geometrical  parameters  of  a  cryocooler.  For  example,  for  an 
ideal  Stirling  cycle  operating  between  T,,  and  Tg,  with  a  volume 
compression  ratio,  r^,  and  using  a  perfect  gas  with  gas 
constant  R,  the  heat  removed  from  the  refrigerated  volume  is 
written  (Huang,  1976:304-305) 

=  R  ln{r^)  (1) 

The  net  work  input  to  the  ideal  cycle  is  written 

=  R  {T,  -  TJ  ln{r,)  (2) 
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To  account  for  the  imperfection  in  the  regenerator,  the 
shortfall  in  heat  transfer  during  one  blow  through  the 
regenerator  must  be  deducted  from  the  ideal  heat  transfer  from 
the  refrigerated  space. 


!2e  =  C’ei  -  (1  -  6)  Qj,  (4) 

where  e  is  the  effectiveness,  i.e.  the  ratio  between  the  heat 
transfer  to  the  regenerator,  Q^,  and  the  ideal  heat  transfer 
to  the  regenerator  under  the  same  conditions .  0^  is  written 


Or  Cy  {T^  Tg) 


(5) 


Hence,  the  COP  for  the  non-ideal  cycle  is 


COP  =  COPj^  -  (1  -  e)  {T^  - 


i?  (Tg  -  Tg)  In(rg) 


(6) 


But,  a  simplified  analysis  shows  (Urielli  and  Berchowitz, 
1984:118) 


e 


N, 


tu 


Ntu  +  1 


(7) 


where  is  the  number  of  thermal  units.  Consequently,  the 
COP  is  written 
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-  T,) 

+  1) 


(8) 


COP  =  COPj^ 


where  A*  =  c^,/ (T^  -  T^) /In  (r^) /R.  But,  for  a  perfect  gas, 
(Tc  -  Tg)  =  Ap  /p/R,  therefore 


COP  -  COP.  -  A**  - - - 

Ap  Ap 


(9) 


where 


A**  =  AVp/R 


From  the  definition  of  the  compactness  factor,  CF  =  f ,  the 

ratio  iV(-u/Ap  can  be  written 


tu  _  D  * 


Ap 


=  B*  CF 


(10) 


where 


and 


B*  =  A^^^/A/Pr2^V(0.5pU-/) 

=  surface  area  of  the  matrix  -  m^ 

A  =  cross  sectional  area  of  the  matrix 
Pr  =  Prandtl  number  (Cp  p/k) 

Ug  =  gas  velocity 


m 


Hence,  the  COP  is  written 


COP  =  COP^  - 


A" 


B*  CF  + 


(11) 


Ap 


This  equation  shows  that  the  best  efficiency  (Carnot)  can  be 
approached  for  large  values  of  CF  and  small  Ap,  although  the 
term  involving  CF  is  usually  much  larger  than  the  inverse  of 
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the  pressure  drop.  The  importance  of  the  heat  transfer  and 
pressure  drop  characteristics  of  a  regenerative  heat  exchanger 
is  explained  more  fully  in  the  next  section. 

Performance  Factors.  There  are  three  factors  which 
significantly  affect  the  COP  of  a  regenerative  refrigeration 
cycle:  1)  regenerator  effectiveness,  2)  pressure  drop,  and  3) 
dead  volume.  All  three  factors  are  influenced  by  the  design 
of  the  regenerator. 

Regenerator  Effectiveness.  Below,  a  definition  for 
regenerator  effectiveness  is  given,  as  well  as  the  importance 
of  effectiveness  for  a  good  COP.  The  most  common  design 
approach  for  maximizing  effectiveness  is  described  and  a 
typical  range  of  effectiveness  values  is  included  to 
facilitate  comprehension  of  current  capabilities. 

The  regenerator  effectiveness  is  defined  as  the  ratio  of 
the  energy  transferred  during  a  blow  through  the  regenerator, 
to  the  energy  that  would  have  been  transferred  in  an  ideal 
regenerator  operating  between  the  same  two  temperatures,  T^ 
and  Tj, .  During  the  cold  blow,  for  example,  the  regenerator 
ideally  receives  gas  at  T^  and  gives  it  up  at  T,,,  meaning  an 
ideal  energy  transfer  per  unit  mass  of  Cp(Tg-T,,),  where  Cp  is 
the  specific  heat  of  the  gas  at  constant  pressure.  However, 
real  regenerators  have  finite  heat  capacities,  and  inducing 
heat  transfer  between  the  gas  and  the  regenerator  matrix 
requires  a  temperature  difference.  Consequently,  the  cold 
blow  gas  leaves  the  regenerator  at  a  temperature  lower  than  T^ 
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and  the  hot  blow  gas  leaves  the  regenerator  at  a  temperature 
higher  than  Tg.  For  high-efficiency,  space-based  cryocoolers, 
a  typical  regenerator  effectiveness  exceeds  0.95. 

The  COP  of  a  Stirling  cycle  is  strongly  dependent  on  the 
regenerator  effectiveness.  In  Fig.  3,  the  results  of  a  study 
by  Tailor  and  Narayankhedkar  (1988)  show  the  COP  plotted 
versus  the  volume  compression  ratio,  which  is  defined  as 
the  ratio  of  the  volume  of  the  gas  before  compression,  to  the 
volume  of  gas  after  compression.  For  a  range  of  compression 
ratios  common  to  cryocoolers,  a  COP  reduction  of  nearly  20% 
occurs  for  just  a  1%  reduction  in  effectiveness.  Clearly,  it 
is  important  to  have  as  high  a  regenerator  effectiveness  as 
possible . 


Figure  3 .  COP  vs  Compression.  Ratio  emd  Effectiveness 

(Tailor  and  Narayankhedkar,  1988:40) 
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High  regenerator  effectiveness  is  achieved  in  the 
following  way.  Small -diameter  wire  (0.02  -  0.20  mm),  closely- 
packed  (100-400  mesh) ,  stainless  steel  or  phosphor-bronze 
screens  are  used.  The  screens  are  punched  into  disks  the 
diameter  of  the  conduit  between  the  expansion  and  compression 
spaces.  A  number  of  randomly-oriented  disks  are  then  stacked 
together  to  form  the  regenerator.  The  tortuous  path  the  gas 
must  travel  through  the  tightly-spaced,  randomly-oriented  wire 
screens  causes  the  gas  to  have  extensive  contact  with  the 
regenerator  material  in  short  axial  distances,  i.e.  the 
regenerator  has  a  large  surf ace-area-to-volume  ratio.  The 
turbulent  flow  has  a  large  heat  transfer  coefficient,  usually 
greater  than  5000  W/m^/°C.  A  linear  temperature  profile 
exists  in  the  regenerator  (Yuan  et  al,  1992;  Tanaka  et  al, 
1990).  Consequently,  as  the  gas  weaves  its  way  through  the 
regenerator,  it  effectively  exchanges  energy  with  the  wire 
screens  which  are  made  from  a  material  with  a  relatively  high 
thermal  capacity.  The  combination  of  large  surf ace-area-to- 
volume  ratio,  large  heat  transfer  coefficient,  and  large 
thermal  capacity  causes  the  regenerator  to  have  a  large 
effectiveness.  Regenerator  configurations  other  than  stacked 
wire  screens  are  used  (Radebaugh  and  Louie,  1986:180; 
Venkatarathnam,  1990),  but  wire  screens  are  the  most  common. 

Pressure  Drop.  A  negative  consequence  of  the 
tightly-packed  porous  material  in  the  regenerator  is  a  large 
pressure  drop,  compared  to  flow  through  an  open  pipe  for 
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example.  Both  inertial  effects  (numerous  accelerations  due  to 
change  in  direction  and  pore  size)  and  surface  shear  stresses, 
contribute  to  this  pressure  drop  (Organ,  1984;  Krazinski, 
1986;  Kays  and  London,  1984;  Armour  and  Cannon,  1968) .  The 
larger  the  pressure  drop  in  the  cryocooler,  the  larger  the  net 
work  input  becomes.  This  reduces  the  COP.  Consequently,  the 
cryocooler  designer  tries  to  minimize  the  losses  associated 
with  pressure  drops  in  the  regenerator. 

However,  minimizing  the  pressure  drop  poses  a  dilemma. 
The  steps  one  would  take  to  reduce  the  pressure  drop  in  the 
regenerator,  e.g.,  making  the  pore  size  large  and 
straightening  the  gas  flow  channels,  would  also  reduce  its 
effectiveness.  As  stated  above,  the  effectiveness  of  the 
regenerator  is  the  dominant  factor  for  determining  cryocooler 
performance  (Martini,  1980:57;  Uriel!  and  Berchowitz ,  1984:99; 
Tanaka  et  al,  1990:177).  Hence,  the  regenerator  designer 
chooses  the  wire  screen  geometry  that  produces  the 
effectiveness  needed  to  fulfill  the  COP  requirement  of  the 
system,  and  the  pressure  drop  plays  a  secondary  role.  But 
there  is  one  other  factor  to  consider. 

The  Dead  Volume.  The  third  factor  which 
significantly  affects  the  performance  of  a  cryocooler  and 
which  is  directly  influenced  by  regenerator  design  is  the  dead 
volume.  The  dead  volume  is  any  gas-filled  space  in  the 
cryocooler  which  is  not  part  of  the  expansion  or  compression 
spaces,  for  example  conduits,  heater  tubes,  and  pore  space  in 
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the  regenerator.  Martini  (1980),  Jones  (1986),  Radebaugh  and 
Daney  (1984) ,  Tailor  and  Narayankhedkar  (1988) ,  and  Romm  and 
Smith  (1993)  have  shown  the  performance  of  the  cryocooler  is 
degraded  by  excessive  dead  volume.  Fig.  4  shows  COP  plotted 
versus  volume  compression  ratio,  r^,  defined  above.  For  a 
typical  of  around  3.0,  and  a  dead  volume  ratio  of  1.0  (here 
defined  as  the  ratio  of  the  dead  volume  to  the  swept  volume  in 
the  expansion  space),  the  COP  is  0.21.  If  the  dead  volume  is 
reduced  by  a  factor  of  25%,  the  COP  can  be  increased  to  0.218, 
which  is  an  increase  of  4%  in  COP.  This  represents  a 
significant  improvement  in  the  performance  of  a  cryocooler 
since  by  the  time  the  compressor  and  solar  panel  efficiencies 
are  included  in  the  overall  energy  savings,  considerable 
weight  and  solar  panel  area  have  been  saved. 

The  porosity  of  the  regenerator  matrix  determines  its 
dead  volume.  Porosity  is  defined  as  the  ratio  of  the  volume 
of  the  dead  space  to  the  total  volume  of  the  regenerator.  The 
value  of  the  porosity  for  most  regenerators  ranges  from  0.5  to 
0.8.  But  these  porosities  also  indicate  over  half  of  the 
space  occupied  by  the  regenerator  is  dead  volume.  If  a  way 
could  be  found  for  achieving  the  effectiveness  needed  for 
system  COP  requirements  without  increasing  the  pressure  drop 
while  also  decreasing  the  dead  volume  in  the  regenerator, 
overall  cryocooler  performance  would  be  improved. 

Rolling  the  wire  screens  would  reduce  their  thickness. 
Stacking  thinner  screens  together  reduces  regenerator  volume 


17 


Volume  compression  ratio, /e 


Figure  4.  Effect  of  Dead  Volume  on  COP 

(Tailor  and  Narayankhedkar ,  1988:41) 


without  reducing  the  volume  of  wire.  Hence,  pore  sizes  in  the 
reduced- thickness  regenerator  are  smaller,  and  the  amount  of 
dead  volume  is  decreased.  The  objective  of  this  research  is 
to  investigate  the  effect  of  rolled  screens  on  the  heat 
transfer  and  pressure  drop  characteristics  of  stacked  wire- 
screen  regenerators  and  to  determine  if  the  reduction  in 
thickness  causes  a  decrease  in  the  compactness  factor. 
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Justification 


It  is  impossible  to  predict  how  the  heat  transfer  and 
friction  characteristics  would  change  if  the  geometry  of  the 
screen  wires  were  flattened.  The  flow  of  gas  through  a  wire- 
screen  regenerator  is  very  complicated.  A  number  of  small- 
scale  phenomena,  such  as  increased  apparent  conductivity, 
channelling  (non-uniform  flow  through  the  porous  medium) ,  and 
secondary  vortices,  influence  the  effectiveness  and  pressure 
drop  of  the  regenerator  (Hutchinson  and  Ross,  1987;  West, 

1986)  .  These  phenomena  are  not  well  understood  nor  quantified 
(Organ,  1992).  If  they  were,  the  precise  geometry  and 
orientation  of  the  regenerator  matrix  could  be  determined  in 
advance.  And  yet,  there  is  no  reason  to  believe  that  screens 
made  from  perfectly  round  wire  give  the  best  results.  They 
are  certainly  the  easiest  to  manufacture,  but  flattening  the 
screens  by  15  or  30  percent,  may  have  only  a  minor  effect  on 
important  flow  characteristics  in  the  regenerator,  like 
pressure  drop  and  heat  transfer  coefficient.  All  the  while, 
the  benefit  of  a  reduction  in  dead  volume  would  be  realized. 

Evidence  exists  that  supports  the  contention  flattening 
the  screens  improves  system-level  performance.  The  first 
piece  of  evidence  is  shown  in  the  top  part  of  Fig.  5.  Cycle 
efficiency  and  indicated  power  for  a  newly-developed  engine 
are  plotted  as  a  function  of  engine  speed  (Nagawa,  et  al, 

1987)  .  These  results  are  for  a  power  generation  cycle,  as 
opposed  to  a  refrigerator,  but  the  conclusions  can  be  directly 
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applied.  The  indicated  efficiency,  which  is  the  ratio  of  the 
indicated  power  output  to  the  net  rate  of  energy  input,  is 
reportedly  improved  by  as  much  as  5%  using  a  regenerator  which 
is  made  of  alternately-stacked,  rolled  (assumed  the  same  as 
pressed)  and  round  wire  screens. 

Other  data  given  in  Fig.  5  indicate  the  pressure  drop 
characteristics  of  a  regenerator  change  if  the  geometry  of  the 
screens  changes.  The  bottom  part  of  Fig.  5  shows  the  pressure 
drop  for  rolled  wire  screens  is  significantly  higher  than  the 
pressure  drop  for  normal  round  wire  screens .  Also  shown  in 
the  bottom  part  of  Fig.  5  is  that  alternately  rolled  and  round 
wire  screens  stacked  randomly  curiously  results  in  a  lower 
pressure  drop.  These  data  indicate  for  the  operating 
conditions  in  the  test,  a  combination  of  rolled  and  round  wire 
screens  reduced  pressure  drop  through  the  regenerator,  while 
using  rolled  wire  screens  alone  increased  it.  However,  a 
range  in  pressure  drops  might  be  achieved  by  rolling  the  wire 
screens  by  various  amounts.  An  optimal  condition  for  overall 
cycle  performance  should  exist. 

Although  the  issue  of  regenerator  effectiveness  is  not 
directly  addressed  in  the  Nagawa  article,  the  following  can  be 
concluded.  If  the  effectiveness  had  been  degraded  excessively 
by  the  arrangement  of  rolled  and  round  wire  screens,  the  cycle 
efficiency  would  not  have  been  improved  as  shown.  In  fact, 
the  authors  reported  their  regenerator  not  only  achieved  a 
reduction  of  flow  friction  loss,  but  also  achieved  effective 
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Figure  5.  NS30A  Performance  and  Pressed  Screen  Results 

(Nagawa  et  al . ,  1987:1799) 
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heat  transfer  (Nagawa  et  al . ,  1987:1799).  Hence,  although  it 
is  not  clear  which  if  any  of  the  three  performance  factors 
discussed  earlier  caused  the  improvement  in  the  engine 
performance,  using  a  mixture  of  rolled  and  unrolled  wire 
screens  in  the  regenerator  altered  the  performance 
significantly.  For  the  NS30A  engine,  the  performance  was 
improved,  but  this  might  not  always  be  the  case. 

The  evidence  given  above  cannot  be  considered  complete. 
Details  about  the  amount  and  type  of  rolling,  the  arrangement 
of  the  screens,  definition  of  terms,  and  the  experimental 
procedure  are  not  available  for  scrutiny.  It  is  not  obvious 
staggering  rolled  and  unrolled  wire  screens  would  cause  an 
improvement  in  pressure  drop.  Perhaps  a  quirk  in  the 
definition  of  terms  or  testing  procedure  caused  the  reported 
results.  Yet,  some  rationale  may  exist  to  explain  these 
improvements  in  efficiency  and  pressure  drop.  Hence,  although 
the  evidence  given  above  seems  to  contradict  the  hypothesis 
given  in  Chapter  I,  it  is  for  only  one  case,  and  the 
documentation  is  not  complete.  In  any  case,  direct 
correlation  between  the  improvement  in  system-level 
performance  and  the  characteristics  of  the  regenerator  as 
measured  by  the  compactness  factor  have  not  been  made.  Since 
a  higher  compactness  factor  correlates  with  higher  engine 
efficiency  (shown  earlier) ,  the  improvement  in  engine 
efficiency  due  to  a  reduction  in  dead  volume  caused  by 
manufacturing  the  regenerator  with  pressed  screens  may 
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outweigh  the  decrease  in  compactness  factor  due  to  the  pressed 
screens  shown  in  this  study.  The  focus  of  the  current 
research  is  on  a  component-level  characteristic,  i.e.,  the 
effect  of  rolled  screens  on  the  compactness  factor,  and  not  on 
the  system- level  performance. 

Other  evidence  exists  which  directly  addresses  the 
hypothesis.  An  experimental  result  by  Wiese  (1993)  using  a 
single-blow,  transient  procedure,  showed  that  the  screens 
could  be  reduced  in  thickness  by  as  much  as  3  0%  without 
significantly  increasing  the  pressure  drop.  The  heat  transfer 
results  from  Wiese's  study  are  not  as  conclusive.  His 
research  indicated  that  at  very  low  Reynolds  number  (based  on 
hydraulic  diameter,  RedH=  Ugdj^/v)  ,  the  heat  transfer 
coefficient  is  significantly  reduced,  while  at  higher  Re^H,  the 
effect  on  heat  transfer  coefficient  is  less  pronounced.  These 
results  indicate  that  some  Reynolds  number  range  may  exist 
where  both  the  pressure  drop  and  heat  transfer  characteristics 
of  the  rolled  screen  regenerator  are  only  changed  slightly, 
e.g.,  the  compactness  factor  is  not  changed,  while  the 
reduction  in  dead  volume  is  significant.  Wiese's  testing 
range  was  for  Re^j^  <  100.  This  range  is  low  compared  to 
typical  cryocooler  operating  conditions  which  is  50  <  Re^jjj  < 
800  (Organ,  1992:  App.  VI;  Seume  and  Simon,  1986:  533-538; 
Walker,  1983:127). 

The  evidence  given  above  seems  to  contradict  the 
contention  that  using  rolled  screens  in  the  regenerator 
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decreases  the  compactness  factor.  Detailed  studies  which 
cover  a  range  of  operating  conditions  and  regenerator 
configurations  are  scarce.  A  literature  search  which 
encompassed  the  last  40  years,  170  documents  (books,  papers, 
contractor  reports) ,  and  research  on  three  continents  (North 
America,  Asia,  and  Europe)  revealed  only  the  Wiese  and  Nagawa 
citations  mentioned  above.  Wiese  described  his  recent 
discussions  with  cryocooler  experts  who  attended  the  8th 
International  Cryocooler  Conference  in  June  of  1994  ,  and  who 
had  attempted  to  improve  regenerative  refrigeration  cycles 
with  rolled  screens,  revealed  mixed  results.  One  researcher 
noted  a  significant  improvement  in  COP  using  rolled  wire 
screens  (as  much  as  30%) ,  while  another  reported  no 
improvement  at  all  (Wiese,  1994) .  A  more  comprehensive  study 
is  needed.  The  study  should  include  varying  the  amount  of 
rolling  (reduction  factor) ,  different  mesh  sizes  and 
porosities,  and  a  range  of  Re  applicable  to  cryocoolers .  A 
better  understanding  of  the  effects  of  screen  thickness 
reduction  on  the  compactness  factor  over  a  wide  range  of 
operating  conditions  would  allow  regenerator  designers  to 
perform  trade  off  studies  between  the  benefits  of  lower  dead 
volume  and  the  detriments  of  lower  compactness  factor. 

In  summary,  regenerative  refrigeration  cycles  are  an 
important  technology.  The  performance  of  these  cycles  is 
improved  by  reducing  the  dead  volume  in  the  system.  Rolling 
the  screens  which  are  used  to  make  the  regenerator  in  these 
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systems  would  reduce  the  dead  volume.  Currently,  the  effects 
of  this  screen  thickness  reduction  on  the  heat  transfer  and 
friction  characteristics  of  the  regenerator,  or  their  ratio, 
the  compactness  factor,  are  neither  quantified  nor  well 
understood.  Hence,  the  hypothesis  given  in  Chapter  I  which 
motivates  the  current  research  was  formulated. 
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Ill .  Methodology 


This  chapter  describes  the  method  for  proving  the 
hypothesis  stated  in  Chapter  I.  The  theory  of  the  classical 
approach  for  finding  the  heat  transfer  coefficient  in  a  porous 
medium  is  given,  and  how  this  approach  was  modified  for  the 
current  research  is  explained.  Next,  the  numerical  model  and 
procedure  for  determining  the  heat  transfer  coefficient  and 
other  important  results  are  described.  The  experimental 
apparatus  and  data  reduction  technique  are  described  in 
Chapter  IV. 

Theory 

This  section  of  the  dissertation  addresses  three  topics. 
First,  a  description  is  given  of  the  step-change  transient 
technique  which  has  been  used  to  determine  heat  transfer 
coefficient  for  porous  material  heat  exchangers  for  nearly 
thirty  years.  Second,  an  explanation  of  the  shortcomings  of 
the  technique  is  presented.  Third,  improvements  to  this 
approach  are  offered  along  with  the  rationale  behind  the 
modifications . 

Step-Change  Transient  Technique.  This  technique  was 
originally  described  by  Pucci,  et  al .  (1967).  As  an 
overview,  the  heat  transfer  coefficient  is  determined  with  the 
step-change  transient  technique  by  comparing  an  experimental 
measurement  of  the  maximum  slope  of  the  outlet  temperature 
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trace  for  a  gas  flowing  through  a  porous  medium  subjected  to 
an  impulsive  change  in  inlet  temperature,  with  an  analytical 
solution  for  the  maximum  slope.  Directly  measuring  the  area, 
energy  flow,  and  temperature  difference  which  define  the  heat 
transfer  coefficient  cannot  be  done  since  instrumenting  the 
tiny  porous  structure  in  a  regenerator  causes  changes  in  the 
local  flow  conditions.  Even  if  a  local  heat  transfer 
coefficient  could  be  determined,  it  would  not  be 
representative  of  the  entire  porous  medium  since  the  velocity 
of  the  fluid  and  the  surface  area  of  the  porous  medium  change 
significantly  from  one  location  to  another.  Hence,  certain 
assumptions  are  made  in  a  simplified  model  of  the  flow,  e.g. 
that  the  heat  transfer  coefficient,  h^onv/  is  constant,  and  the 
magnitude  of  the  heat  transfer  coefficient  is  determined  by 
comparing  the  results  of  an  easily  measured  quantity  like  the 
outlet  temperature,  to  the  analytical  result. 

Predicting  the  magnitude  and  direction  for  flow  of  a  gas 
through  a  porous  medium  regenerator  is  complicated  by  the 
numerous  accelerations  and  changes  in  direction  experienced  by 
the  gas  as  it  winds  its  way  through  the  randomly  oriented  pore 
structure.  To  make  the  task  manageable,  certain  assumptions 
are  made  about  the  flow.  The  pore  sizes  are  assumed  uniform 
and  the  flow  is  assumed  to  be  one-dimensional.  The  average 
velocity  at  any  location  in  the  regenerator  is  presumed  to  be 
in  the  direction  of  the  regenerator  axis  with  a  magnitude 
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equal  to  the  total  mass  flow  rate  divided  by  an  average 
density,  divided  by  the  average  cross-sectional  area  for  flow. 
If  the  Mach  number  is  kept  below  0.2,  compressibility  effects 
can  be  ignored  (Walker,  1983: Part  I,  Section  5.3)  and  the  mass 
velocity,  Ug,  is  constant.  If  it  is  further  assumed  that 
temperature  and  pressure  differences  throughout  the 
regenerator  are  kept  small  (say  <  20°  C  and  50  kPa) ,  then  gas 
properties  are  constant.  Further,  if  convective  heat  transfer 
between  the  gas  and  the  regenerator  is  the  only  energy 
exchange  mechanism  considered,  the  energy  equations  for  the 
gas  and  matrix  are  written  as  follows: 

Gas : 

dT  dT 

-  ^  Cp  Ax  +  hcoNv^s  =  p  A  $  Ax  (12) 


Matrix: 


^CONV 


Pw  (1  - 


dT^ 

dt 


(13) 


where  iti 

Cp,  Cy 


Cpi 

Ax 

^conv 


is  the  mass  flow  rate  (kg/sec) 
is  the  specific  heat  of  the  gas  at 
constant  pressure  and  volume, 
respectively  (J/kg/°C) 
is  the  specific  heat  of  the  matrix 
(J/kg/°C) 

are  gas  and  matrix  temperature, 
respectively  (°C) 
is  a  linear  increment  (m) 
is  the  heat  transfer  coefficient 
(W/mV°C) 


28 


and 


P/Pm  the  gas  and  matrix  densities 

(kg/m^) 

A  is  the  tube  cross  sectional  area  (m^’ 

is  matrix  surface  area  (m^) 

X,  t  are  the  axial  location  and  time 

O  is  the  porosity  defined  as  the  volume 

of  the  gas  divided  by  the  total 
volume  of  the  regenerator. 


Pucci,  et  al .  (1967:30)  showed  that  with  the  introduction  of 
the  appropriate  dimensionless  variables,  along  with  an 
adiabatic  boundary  condition  and  the  impulsive  step -change 
initial  condition,  this  equation  could  be  solved  for  the 
temperature  of  the  gas  as  a  function  of  time  and  axial 
location.  The  solution  is  an  infinite  series  of  Bessel 
functions  which  depends  on  a  previous  knowledge  of  the  heat 
transfer  coefficient.  However,  Locke  (1950)  showed  that  the 
maximum  slope  of  the  dimensionless  temperature  at  the  outlet 
is  solely  a  function  of  the  heat  transfer  coefficient: 


d&  I 
Imax 


NtU^ 
s/NtU  T 


[  -i  (2  i  ^Ntu  t)  ]  exp  ■  ^ 


(14) 


where 

and 


0  is  a  dimensionless  temperature 

X  is  a  dimensionless  time 

is  the  number  of  thermal  units  = 

{h^onv  As/m/Cp) 


This  solution  to  the  governing  equations  of  the  step-change 
transient  approach  offered  an  elegant  method  for  determining 
the  heat  transfer  coefficient.  First,  introduce  an  impulsive 
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step-change  in  inlet  temperature  to  the  flow  of  gas  through  a 
porous  medium  regenerator.  Next,  measure  the  maximum  slope  of 
the  outlet  temperature  trace.  Finally,  from  the  maximum 
gradient  of  the  dimensionless  outlet  temperature  0  with 
respect  to  the  dimensionless  time,  a  value  for  Ntu  can  be 
calculated.  From  this  Ntu,  the  hcow  can  be  calculated.  This 
method  for  determining  the  heat  transfer  coefficient  will  be 
referred  to  as  the  analytical  method  due  to  Pucci,  et  al . ,  in 
the  remainder  of  this  text. 

Other  researchers  take  a  slightly  different  approach. 
Liang  and  Yang  (1975)  used  outlet  temperature  data  at  four  or 
five  points  to  calculate  a  minimum  root-mean-squared 
difference  between  an  experimental  outlet  temperature  and  one 
determined  by  an  analytical  solution.  In  fact,  since  its 
introduction,  several  modifications  and  improvements  have  been 
made  to  the  step- change  transient  technique  (Stang  and  Bush, 
1974;  Yagi,  1991),  but  the  basic  approach  remains  the  same. 

Problems  with  the  Step-Change  Transient  Technique.  There 
are  three  major  problems  with  the  technique  described  in  the 
last  section.  First,  obtaining  an  adequate  approximation  to 
an  exact  step-change  in  the  laboratory  is  impractical. 
Second,  if  one  chooses  to  use  the  maximum  slope  approach,  it 
is  very  difficult  to  obtain  an  accurate  slope  from 
experimental  data.  Third,  the  model  ignores  some  important 
physics.  In  this  section,  the  cause  of  these  discrepancies 
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between  the  theory  and  experimental  results  is  reported  as 
well  as  the  effect  on  the  determination  of  the  heat  transfer 
coefficient . 

With  regard  to  the  first  problem,  any  attempt  at 
producing  a  step-change  in  temperature  in  the  laboratory  will 
have  some  time  delay.  Consequently,  whether  using  a  maximum 
slope,  or  a  root-mean-squared  difference,  the  calculated  heat 
transfer  coefficient  will  be  smaller  than  it  should  be, 
because  the  output  temperature  trace  will  be  less  steep  than 
it  should  be  if  the  input  temperature  trace  is  less  steep  than 
ideal . 

The  second  major  problem  is  associated  with  using  the 
maximum  slope.  Digitally  acquired  experimental  data 
inherently  lack  smoothness.  Digital  data  are  not  continuous, 
a  mathematical  requirement  for  the  existence  of  the  slope. 
Data  are  collected  at  sampling  rate  intervals  which  are  small, 
but  are  not  continuous,  even  in  the  limit.  Second, 
experimental  data  have  some  scatter.  Gradients  can  shift 
wildly  from  one  set  of  sampling  points  to  the  next, 
particularly  if  one  uses  a  very  high  sampling  rate  and  the 
time  step  between  data  points  is  small.  One  could  use 
filtering  techniques  or  Fourier  transforms  to  improve  the 
appearance  of  the  data,  but  these  techniques  have  a 
deleterious  effect  on  the  final  value  for  the  slope.  For 
example,  large  order  filters  smooth  out  any  discontinuities  in 
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the  data,  but  this  can  significantly  change  the  value  of  the 
slope.  Other  researchers  have  noted  these  problems  also 
(Shoup,  1977:  220). 

The  third  major  problem  with  the  step-change  technique  is 
that  heat  transfer  to  the  tube  wall  is  ignored.  The 
importance  of  the  wall  became  apparent  during  check-out  of  the 
experimental  set-up.  During  some  preliminary  transients,  an 
unexpected  outlet  temperature  behavior  was  observed.  In 
Fig.  6,  a  temperature  transient  for  an  empty  tube  is  shown. 
The  inlet  and  outlet  temperatures  at  the  start  of  the  test  run 
were  at  the  same  steady-state  value.  When  the  transient 
began,  the  inlet  temperature,  T,^,  abruptly  fell,  and  soon 
began  to  relax  towards  another  steady-state  value.  Similarly, 
the  outlet  temperature  trace,  Tg2,  fell  sharply,  and  relaxed 
toward  a  new  steady-state  value.  But  the  later  time  Tg2  is 
approximately  1-2°  C  above  that  of  the  later  time  This 

indicated  a  discrepancy  in  the  energy  exchange  terms  included 
in  Pucci's  model  in  Eqs  (12)  and  (13).  and  should 

relax  toward  the  same  late-time  temperature  since  only 
interaction  between  the  gas  and  the  matrix  is  relevant. 

Before  taking  any  heat  transfer  coefficient  data,  the 
discrepancy  had  to  be  explained.  After  repeated  careful 
calibrations  of  the  thermistors  which  measure  Tg2_  and  and 

modifications  to  the  instrumentation,  the  temperature  gap  at 
late  time  remained.  Apparently,  the  tube  and  its  connections 
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Traces  Showing  Temperatures  for  Empty  Tube 


Figure  6.  Temperature  Gap  for  Empty  Tube 


were  supplying  enough  energy  to  measurably  influence  the 
temperature  of  the  gas  for  a  significant  time  after  the 
transient  was  completed.  A  simple  calculation  shows  the 
thermal  inertia  of  the  tube  and  environment  is  much  larger 
than  the  thermal  inertia  of  the  matrix.  If  the  tube  and 
matrix  are  treated  in  a  lumped-parameter  sense,  and  if  the 
heat  transfer  coefficient  between  the  gas  and  the  tube  is 
considered  the  same  as  between  the  gas  and  the  matrix,  the 
ratio  of  time  constants  for  heat  transfer  for  a  length  of  tube 
to  that  for  a  length  of  matrix  is  written: 
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(1  -  <E>) 


(15) 


where  Tj.  is  the  time  constant  for  a  length  of  tube 

(sec) 

Xr  is  the  time  constant  for  a  length  of 
regenerator  matrix  (sec) 
tf.  is  the  tube  thickness  (m) 

O  is  the  porosity 

and  is  the  specific  surface  area  for  heat 

transfer  per  unit  volume  of  the 
matrix  (m‘^) 

For  the  twelve  SUs  used  in  the  current  research,  this  ratio 
ranged  from  16.4  for  SU  #1  (100  mesh.  Reduction  Factor  =  1.0), 
to  47.2  for  SU  #12  (250  mesh,  RF=0.5)  .  Consequently,  the  heat 
transfer  in  the  matrix  occurred  at  a  higher  rate  than  that  for 
the  tube.  Thus,  the  early-time  heat  transfer  should  be 
studied  to  get  the  best  estimate  of  the  heat  transfer 
coefficient  for  the  matrix.  However,  the  step-change 
transient  technique  concentrates  on  relatively  late  time 
events.  Hence,  both  the  maximum  slope  and  the  minimum  root- 
mean-squared  difference  approach  to  the  classical  step-change 
transient  technique  result  in  erroneous,  low  values  for  the 
heat  transfer  coefficient. 

In  summary,  the  step-change  transient  technique  is  a 
simple  approach  for  determining  the  heat  transfer  coefficient 
for  a  complicated  phenomena,  flow  through  a  porous  media. 
However,  it  focuses  on  properties  of  the  outlet  temperature 
trace  to  calculate  the  heat  transfer  coefficient.  If  these 


34 


properties  are  altered  by  the  inability  to  accurately  produce 
a  step-change  in  temperature,  or  by  the  measurement  of  the 
maximum  slope,  or  by  heat  transfer  between  the  gas  and  the 
tube  wall,  then  an  erroneous  value  for  the  heat  transfer 
coefficient  will  be  calculated.  An  improved  approach  is 
offered  next,  based  on  the  same  fundamentals  of  the  step- 
change  technique,  but  which  mitigates  the  problems  discussed 
above . 

Improvements  to  the  Step-Change  Transient  Technique. 
Three  major  modifications  to  the  classical  method  described 
above  are  made  in  the  current  research.  The  first  is  to 
eliminate  the  assumption  that  an  exact  step-change  in 
temperature  occurs  at  the  inlet.  This  is  accomplished  by 
measuring  the  inlet  temperature  trace  and  using  it  in  a 
numerical  model.  Current  data  acquisition  systems  and 
computers  make  this  task  manageable. 

The  second  major  improvement  is  the  inclusion  of  the  tube 
which  surrounds  the  matrix  in  the  model.  Before  deciding  on 
just  how  to  do  this,  further  investigation  on  the  effects  of 
the  tube  was  done.  SU  #1  was  instrumented  with  a  small 
thermocouple  on  the  outside  surface  of  the  tube  midway  between 
the  ends .  The  temperature  of  the  tube  was  monitored  during  a 
transient.  The  results  are  shown  in  Fig.  7.  The  tube 
temperature  stayed  constant  at  its  original  value  until  midway 
through  the  transient,  then  it  rose  rapidly  and  relaxed  to  a 
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steady  state  value  below  that  of  either  the  inlet  or  outlet 
gas . 

Three  comments  about  transient  heat  transfer  to  the  tube 
are  important.  First,  the  data  show,  when  a  regenerator  is  in 
the  flow,  a  temperature  wave  travels  down  the  tube  and  matrix, 
i.e.  rather  than  a  smooth  transition  between  the  initial  and 
final  step-change  in  temperature,  the  transition  is  abrupt, 
and  the  location  of  this  abrupt  transition  moves  down  the  axis 
of  the  regenerator  with  time.  Second,  because  of  this 
traveling  temperature  wave,  whether  the  tube  is  a  thin-walled 
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conductor  or  an  insulator,  the  boundary  condition  for  the  tube 
changes  with  time.  No  simple  analytic  boundary  condition, 
e.g.  adiabatic  or  isothermal,  could  replace  the  effect  of  the 
tube.  Third,  a  comparison  between  the  time  constant  for  a 
unit  length  of  tube  to  that  for  a  unit  length  of  matrix  was 
shown  earlier.  This  comparison  showed  that  the  tube  by  itself 
has  a  much  longer  relaxation  time  than  the  matrix.  The  tube 
is  connected  on  either  end  to  copper  tubing.  The  effect  of 
the  adjoining  web  of  highly  conductive  copper  tubing  is  to 
increase  the  relaxation  time  of  the  tube  temperature  by 
conducting  heat  through  the  sealed  connections  at  each  end. 
This  effect  of  conduction  through  the  ends  of  the  tube  would 
require  additional  instrumentation  and  modelling.  In  view  of 
the  above  considerations,  a  simple  lumped-parameter  approach 
was  assumed  for  inclusion  of  the  tube,  i.e.  the  heat  transfer 
coefficient  between  the  gas  and  the  tube  is  assumed  to  be  the 
same  as  between  the  gas  and  the  matrix,  and  the  tube  wall 
temperature  is  assumed  constant  in  the  radial  and  azimuthal 
directions  at  each  axial  location  of  the  tube.  Since  the 
choice  of  heat  transfer  coefficient  depends  on  early  time 
criteria  (see  below) ,  and  the  time  constant  for  the  matrix  is 
at  least  16-40  times  faster  than  the  time  constant  for  the 
tube  and  surroundings,  the  results  obtained  by  treating  the 
tube  as  a  lumped  thermal  mass  are  acceptable.  The  merits  of 
this  approach  are  discussed  below  and  in  Chapter  IV. 
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The  third  major  improvement  deals  with  the  criteria  for 
selecting  the  heat  transfer  coefficient.  A  characteristic  of 
porous  media  is  that  it  acts  like  a  thermal  sponge.  The  rest 
of  this  section  describes  this  characteristic  in  more  detail, 
defines  a  parameter  of  each  SU  called  the  sponge  effect  delay- 
time,  SEDT,  and  explains  how  this  parameter  can  be  used  to 
calculate  the  heat  transfer  coefficient. 

The  sponge  effect  of  a  regenerator  can  be  demonstrated 
from  results  shown  in  Fig.  6  and  Fig.  8.  For  the  empty  tube 
in  Fig.  6,  the  outlet  temperature  trace  changes  almost 
immediately,  in  tandem  with  the  inlet  temperature  trace.  The 
difference  between  the  two  traces  is  attributed  to  some  heat 
transfer  between  the  tube  and  the  gas.  However,  when  a 
regenerator  is  present  as  in  Fig.  8,  the  gas  at  the  outlet 
stays  at  nearly  the  initial  steady  state  temperature  for  some 
time  before  it  begins  to  change.  As  discussed  in  Chapter  II, 
porous  medium  regenerators  are  effective  in  capturing  a  large 
amount  of  energy  from  the  incoming  stream  of  gas .  This 
effectiveness  of  the  regenerator  matrix  causes  the  observed 
delay  time  between  the  beginning  of  the  transient  and  when  the 
outlet  temperature  begins  to  fall. 

The  definition  of  the  sponge  effect  delay  time  is  as 
follows:  It  is  the  time  it  takes  for  the  outlet  gas 
temperature  to  change  a  measurable  amount.  For  the  current 
research,  the  definition  of  a  measurable  temperature  change  is 
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that  the  value  of  the  outlet  temperature  deviates  from  its 
running  average  by  more  than  twice  the  accuracy  of  the 
temperature  measuring  device,  or  about  0.4°  C  for  the  current 
research . 

Using  the  sponge  effect  delay  time  in  conjunction  with 
the  step-change  transient  technique  gives  a  better  value  for 
the  heat  transfer  coefficient  than  either  the  maximum  slope  or 
the  minimum  root-mean- squared  difference  approach  described 
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above.  Three  arguments  can  be  made  to  substantiate  this 
claim. 

First,  the  best  value  for  the  heat  transfer  coefficient 
will  be  found  during  that  time  when  heat  transfer  between  the 
gas  and  the  matrix  dominates.  Early  in  the  transient,  the 
interaction  between  the  tube  and  the  gas  is  less  important 
than  energy  exchange  between  the  gas  and  matrix  because  the 
surface  area  per  unit  length  of  the  porous  medium  ranges  from 
40.5  times  larger  (100  mesh)  to  101  times  larger  (250  mesh) 
than  the  surface  area  of  the  tube  (Miyabe,  et  al . ,  1982:1840) 
while  the  temperature  gradient  is  approximately  the  same.  The 
analysis  of  time  constants  for  a  length  of  tube  given  above 
showed  that  the  tube  and  its  environment  are  sixteen  to  forty 
seven  times  slower  to  respond  to  the  transient  than  the 
regenerator.  Hence,  during  the  sponge  effect  delay  time,  heat 
transfer  between  the  gas  and  the  matrix  dominates  energy 
exchange.  The  assumptions  made  in  the  numerical  model  more 
nearly  match  the  physics  of  the  flow  during  this  early  time, 
and  a  better  result  for  the  heat  transfer  coefficient  is 
obtained . 

The  second  part  of  the  argument  that  a  heat  transfer 
coefficient  based  on  the  sponge-effect  delay  time  would  be 
better  than  one  based  on  either  the  maximum  slope  or  a  minimum 
root-mean- squared  difference  is  that  the  delay  time  is  a 
characteristic  of  the  regenerator  while  the  outlet  temperature 
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trace  is  not.  Fig.  9  shows  calculated  outlet  temperature 
traces  from  the  numerical  model  described  below  for  an 
artificial  step-change  in  inlet  temperature.  One  of  the 
traces  includes  the  effects  of  the  tube,  the  other  shows  the 
result  without  tube  effects.  Clearly,  one  can  obtain  a  wide 
range  of  slopes  for  the  same  regenerator  since  the  temperature 
trace  with  no  tube  effect  is  much  steeper  than  the  trace 
including  the  tube  effect.  The  minimum  root-mean-squared 
difference  approach  is  no  better  since  the  heat  transfer 
coefficient  chosen  with  this  approach  gives  the  best  match 
between  the  experimental  outlet  temperatures  and  those 
calculated  by  a  selected  analytical  model,  whether  the 
analytical  model  includes  all  the  important  physics  of  the 
flow  or  not.  Hence,  an  approach  which  includes  choosing  the 
heat  transfer  coefficient  based  on  the  outlet  temperature 
trace  has  inherent  difficulties  which  are  better  addressed  by 
using  the  sponge  effect  delay  time. 

The  third  part  of  the  argument  for  the  improvement 
offered  by  the  current  method  can  also  be  seen  from  Fig.  9. 
Whereas  a  wide  range  of  outlet  temperature  traces  can  be 
manufactured  depending  on  how  much  the  tube  and  other  factors 
affect  the  flow,  the  delay  time  is  nearly  insensitive  to  these 
extraneous  factors.  For  the  example  in  Fig.  9,  the  maximum 
slopes  differ  for  the  two  cases  by  approximately  150%,  whereas 
the  delay  times  are  only  different  by  about  5%  which  can  be 
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Outlet  Traces  with  and  without  Influence  of  Tube 


Figrure  9.  Comparison  of  Trace  for  Tube  and  No  Tube 


seen  by  comparing  the  time  it  takes  for  the  two  traces  to 
change  by  0.4°C  from  the  initial  temperature.  Hence,  the 
sponge  effect  delay  time  is  an  easily  measured  parameter  of  a 
SU  which  gives  better  results  for  the  heat  transfer 
coefficient  than  other  parameters  described  above. 

To  summarize,  one  of  the  objectives  of  this  research  is 
to  measure  the  heat  transfer  coefficient  and  friction  factors 
for  flow  through  regenerators  made  of  stacked  screens  of 
varying  reduction  factors.  The  heat  transfer  coefficient 
cannot  be  measured  directly,  so  an  experimental  set-up  was 
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built  to  measure  the  inlet  and  outlet  gas  temperature  during 
a  quasi-step-change  transient,  along  with  other  important 
measurements  like  pressure  and  mass  flow  rate.  This  data  was 
put  into  a  form  suitable  for  use  with  a  numerical  model  of  the 
flow.  This  model  generates  an  outlet  temperature  trace  from 
which  a  characteristic  parameter  of  the  flow,  the  sponge 
effect  delay  time,  can  be  determined.  The  criteria  for 
choosing  the  heat  transfer  coefficient  was  the  best  match 
between  the  sponge  effect  delay  time  measured  from 
experimental  Tg2  data  and  the  one  calculated  by  the  numerical 
model.  Details  about  the  experimental  set-up  and  data 
reduction  technique  are  given  in  Chapter  IV.  A  description  of 
the  numerical  model  and  computational  process  is  next. 

Numerical  Model 

In  this  section,  the  numerical  model  for  flow  through  a 
regenerator  is  described,  and  the  process  for  choosing  the 
best  heat  transfer  coefficient  is  given.  The  discussion  is  in 
two  parts.  First,  the  governing  equations,  assumptions, 
boundary  and  initial  conditions,  and  definitions  for  the 
FORTRAN  program  called  tsie  (for  temperature  solution, 
incompressible,  explicit)  are  given.  Second,  the  way  that 
tsie  was  used  within  a  larger  FORTRAN  program  called  main  to 
determine  the  heat  transfer  coefficient  is  described. 

Numerical  Model  of  the  Transient.  The  FORTRAN  program 
named  tsie  is  an  explicit,  finite  difference  numerical  model 
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of  an  quasi-step-change  transient.  Its  model  is  the  same  as 
the  one  used  in  the  analytical  approach  of  Pucci,  et  al .  , 
described  above  for  the  step-change  transient  solution  except 
for  two  things:  1)  the  inlet  temperature  trace  is  used 

rather  than  an  assumed  step-change,  and  2)  the  influence  of 
the  tube  is  included  in  the  model.  A  lumped  parameter 
approach  to  the  gas /tube  interaction  is  used  and  the  heat 
transfer  coefficient  between  the  tube  and  the  gas  is  assumed 
to  be  the  same  as  the  heat  transfer  coefficient  between  the 
gas  and  the  regenerator  matrix.  Other  assumptions  in  the 
model  include  no  axial  conduction  (Wiese,  1993:2-7),  no 
conduction  between  the  tube  and  the  matrix  (they  are  locally 
at  nearly  the  same  temperature),  perfect  gas  relations, 
constant  properties,  low  Mach  numbers,  one-dimensional, 
incompressible  flow,  and  a  constant  heat  transfer  coefficient. 
With  these  assumptions,  the  momentum  equation  shows  the  mass 
flow  rate,  and  consequently,  the  mass  velocity,  is  a  constant. 
The  important  terms  for  the  energy  equation  are  shown  in 
Fig.  10.  There  is  internal  energy  flowing  with  the  gas, 
convective  heat  transfer  between  the  gas  and  the  matrix, 
convective  heat  transfer  between  the  tube  and  the  matrix, 
conduction  in  the  tube,  and  change  in  internal  energy  of  the 
tube,  matrix,  and  gas  in  the  control  volume.  The  outside 
surface  and  ends  of  the  tube  were  considered  adiabatic.  With 
these  interactions  and  assumptions  in  mind,  three  governing 
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TERMS  FOR  THE  ENERGY 


CONDUCTION  IN: 
Ks  Atat  62  (Tt)/8x2 


CONTROL  VOLUME 


TUBE  WALL  1 


HEATTRANSFER,  GASTOTUBE:  hoonv  Ast  (Tg  -  Tt) 


GAS  FLOW  IN: 
•  m  cp  6(Tg)/5x 


HEATTRANSFER,  GAS  TO  MATRIX:  hoonvAs  {Tg  -  Tr) 


OUTSIDE  OF  TUBE  INSULATED 


CONDUCTION  OUT 


GAS  FLOW  OUT 


1  CHANGE  IN  THERMAL  CAPACITANCE  OF  THE  TUBE;  pmAcxt  Axcm  6CTt)/6t 

2  CHANGEN  IN  THERMAL  CAPACITANCE  OF  THE  MATRIX;  pm  A  Ax  (1.0- cm  6(Tr)/6t 

T  CHANGE  IN  THERMAL  CAPACITANCE  OF  THE  GAS:  p  A  Ax  ®  cp  6(Tg)/6t 


Figure  10.  Control  Volxme  for  Energy  Equation 


equations  can  be  written  as  follows: 


Cp  +  -^COW  ^sr  )  ^CONV  ^ST  (  ) 

9  Cy  A  Lx 


(16) 


Regenerator : 


^coNv  ^g  )  Pm  Cm  L  (1  $ )  Ax 
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Tube : 


^cx 


dx^ 


^CONV  ^ST  ( 


T.-T„ 


)  =  P 


M 


^CX 


dt 


(18) 


where 


and 


m 

^p/ 


h 


conv 


P 

A 


O 

Pm 

Cm 

Ax,  At 


is  the  mass  flow  rate  (kg/sec) 
are  the  specific  heat  of  the  gas  at 
constant  pressure/volume , 
respectively  (J/kg/°C) 
is  the  heat  transfer  coefficient 
(W/mV°C) 

is  the  regenerator  surface  area  (m^) 
is  the  tube  surface  area  (m^) 
is  the  cross-sectional  area 
associated  with  the  tube  wall 
thickness  (m^) 

are  the  gas,  regenerator,  and  tube 
temperatures,  respectively  (°C) 
is  the  gas  density  (kg/m^) 
is  the  cross-sectional  area  of  the 
inner  diameter  of  the  tube  (m^) 
is  the  porosity 

is  the  density  of  the  steel  (kg/m^) 
is  the  specific  heat  of  steel 
(J/kg/oC) 

are  the  axial  location  and  time 
increments,  respectively  (m, sec) 


By  grouping  geometrical  and  physical  properties  together, 
using  a  first-order  accurate  forward  difference  in  time,  and 
a  second-order  accurate  central  difference  in  space,  the 
equations  can  be  recast  as  follows: 

Gas  : 


TG: 


'N 


{- 


TG, 


i-i 

'N+l 


TG^-l 


Kj,  +  Kj  +  Kj,  TGm'^)  / 


(19) 
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Regenerator : 


TR^  =  TG^~^  +(1.0  -  K^) 


(20) 


Tube : 

rr/  =  (TT^ll  +  TT^:l  )  +  Kp  Tr^"" 


where  I,N  is  the  time, space  index 

=  2.0  St  (As/Af) 

Ks  =  (2.0/Y)  {C/Ug) 

Kc  - 

Kj  =  2.0  St  {Asr/Af) 

Kj,  -  K^-K^-Kj 

Ke  —  CLe  a  t 

Kp  =  At/Xt, 

Ke  =  1.0- (2.0  Ke)-Kp 

Ap  =  ^  A  (flow  area) 

St  =  Nu/Re/Pr  (Stanton  number) 

C  =  Courant  Number,  Ax/At 

ttt  =  thermal  diffusivity  of  the  tube 

Y  —  Cp/ Cy 

are  thermal  time  constants  for  the 
matrix  and  tube 

and  Ug  is  gas  mass  velocity 

With  the  finite  difference  form  of  the  governing 
equations  given,  the  next  step  was  to  define  the  initial  and 
boundary  conditions.  The  initial  condition  is  the  same  for 
each  of  the  three  equations  above.  The  gas,  regenerator,  and 
tube  are  all  at  the  same  initial,  steady-state  temperature. 
This  temperature  was  determined  by  an  averaging  of  the 
experimental  gas  temperatures  for  the  one  hundred  time  steps 
before  the  transient  begins.  The  program  begins  calculations 
at  t=Xo  and  continues  calculating  until  t=Tf,  the  end  of  the 
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transient.  There  are  two  boundary  conditions:  1)  the  inlet 
gas  temperature  equals  the  experimentally  measured  value  of 
and  2)  there  is  no  conduction  from  the  ends  of  the  tube, 
i.e.  an  adiabatic  boundary  condition. 

Calculations .  A  printout  of  the  computer  code  is 
given  in  Appendix  C.  A  summary  of  the  way  tsie  calculates  the 
outlet  temperature  trace,  T,2,  is  given  here . 

Step  #1:  Certain  variables  are  input  to  the  model,  i.e. 
the  starting  guess  for  the  heat  transfer  coefficient,  h^onv 
the  maximum  number  of  space  increments,  NMAX,  and  the  value  of 
the  time  increment,  At. 

Step  #2  :  The  geometrical,  parametrical,  and  experimental 
data  for  a  given  test  run  are  read  into  the  program  from  the 
data  file  for  that  test  run. 

Step  #3:  The  initial  temperature,  Tq,  is  calculated. 

The  gas,  tube,  and  regenerator  are  set  equal  to  Tq  at  the 
initial  time  step,  TO. 

Step  #4:  The  values  of  the  coefficients  of  the  finite 
difference  equation  are  calculated. 

Step  #5:  The  temperatures  for  the  current  time  step  at 
each  spacial  location  of  N  are  calculated  with  the  finite 
difference  equations,  Eqs  (19),  (20),  and  (21) .  This  is  done 

using  temperature  values  at  the  current  (I)  and  one  previous 
time  step  [I-l] ,  and  at  locations  one  space  increment  before 
{N-1)  and  one  beyond  {N+1)  the  current  space  increment. 
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step  #6:  The  outlet  temperatures,  are  saved  to  a 

file. 

Step  #7 :  The  current  time  step  temperatures  are  renamed 
as  previous  time  step  temperatures. 

Step  #8:  The  process  is  repeated  until  the  end  of  the 
transient,  if. 

Stability.  The  FORTRAN  program  tsie  becomes 
unstable  at  certain  combinations  of  Ax  and  At.  This  is  common 
for  explicit  finite  difference  models.  Because  three 
interacting  equations  are  involved  in  the  model,  determining 
the  stable  time  step  involves  finding  the  eigenvalues  for  a 
matrix  of  coefficients  and  setting  the  largest  eigenvalue  less 
than  1.0  for  stability.  This  matrix  would  be  large  (NMAX  by 
IMAX)  ,  sparse,  and  different  for  each  value  of  the  heat 
transfer  coefficient  used.  After  reviewing  the  requirements 
for  main  (that  a  high  and  low  guess  for  the  heat  transfer 
coefficient  be  input,  see  below) ,  a  less  complicated  method 
was  sought  to  fulfill  the  stability  requirement.  A 
conservative  approach  is  merely  to  ensure  the  coefficients  in 
the  governing  equations  are  all  positive  (Kreith,  1973:186; 
Patankar,  1980:37) .  For  a  given  value  of  the  ratio  between 
the  time  and  the  space  increment,  or  the  Courant  number,  C, 
the  maximum  heat  transfer  coefficient  for  stability  is 


_  C  Cp  p  Ap 
{A^  +  Y 


(22) 
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where 


and 


Cp  is  the  specific  heat  of  the  gas  at 
constant  pressure  (J/kg/°  C) 

C  is  the  Courant  number,  Ax/At 

p  is  the  gas  density  (kg/m^) 

Ap  is  the  flow  area  in  the  matrix,  <E>  A  (m^) 

Ag  is  the  surface  area  of  the  matrix 

Asj.  is  the  surface  area  of  the  tube  (m^) 

Y  is  the  ratio  Cp/Cy  for  the  gas . 


This  stability  requirement  gives  a  rough  estimate  for  the 
maximum  heat  transfer  coefficient  allowed  for  a  certain  Ax  and 
At.  With  this  criteria,  any  value  for  h^onv  could  have  been 
used,  and  no  stability  problems  were  encountered. 

Validation  of  the  Numerical  Model.  The  model  was 
validated  in  two  ways.  Both  involve  a  comparison  between  a 
known  analytical  solution  for  the  temperature,  and  one 
calculated  by  the  current  numerical  model  (with  some 
modifications  explained  below) . 

In  the  Theory  section  of  this  chapter,  the  exact 
analytical  solution  to  the  classical  step-change  transient 
method  is  presented.  The  current  finite  difference  model 
should  give  the  same  solution  as  the  analytical  one  when  two 
modifications  are  made.  First,  an  artificial,  impulsive  step- 
change  in  temperature  was  fabricated  numerically  to  use  as  the 
input  trace,  Second,  the  area  of  the  tube  surface,  A^, 
was  set  equal  to  zero.  This  effectively  disengages  the 
interaction  between  the  tube  and  the  flowing  gas.  Hence,  all 
the  details  of  the  model  by  Pucci,  et  al . ,  (1967)  match  the 
current  modified  model.  The  results  for  the  maximum  slope  of 
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should  be  the  same  for  both. 


The 


the  outlet  trace,  T,2, 
smooth,  artificially  fabricated  inlet  data  produced  a  smooth 
output  for  which  the  slope  was  determined  by  a  second-order 
accurate,  five  point,  finite  difference  definition  (Anderson, 
et  al . ,  1984:45)  .  This  procedure  was  followed  for  a  test  case 
with  an  arbitrarily  chosen  heat  transfer  coefficient,  an 
artificial,  impulsive  step-change  from  10  to  30°  C,  and  the 
properties  of  SU  #09.  To  determine  the  influence  of  At  and  Ax¬ 
on  accuracy,  the  time  increment.  At,  was  varied  from  10'^  to 
10“’^  seconds,  and  the  number  of  space  increments  was  varied 
from  twenty  one  to  one  hundred  twenty  one.  By  fitting  the 
difference  in  maximum  slope  between  the  analytical  and 
numerical  solutions  for  the  outlet  temperature  to  a  curve,  and 
taking  a  limit  as  NMAX  became  very  large,  the  error  in  maximum 
slope  approached  0.01%.  These  results  are  shown  in  Fig.  11. 
For  the  current  research,  NMAX  of  101  was  deemed  the  best 
balance  between  accuracy  (  1.8%),  and  the  length  of  time  it 
took  to  run  the  numerical  program. 

The  second  test  for  validating  the  numerical  model  also 
involved  using  the  artificial,  impulsive,  twenty  degree  step- 
change  in  inlet  temperature.  In  this  case,  the  temperature 
trace  of  the  matrix  at  the  second  special  increment  was 
observed.  Since  the  wire  diameter  of  the  mesh,  d„,  is  so 
small,  and  the  thermal  conductivity  of  the  304  stainless 
steel,  k^,  is  so  large,  treating  the  matrix  as  a  lumped 
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Figure  11.  Verification  of  Numerical  Approach  for  Large 

NNAX 


thermal  mass  is  permissible.  The  average  Biot  Number  {h^onv 
dw/kg)  for  the  thirteen  SUs  equals  0.048  which  meets  the 
requirement  that  it  be  less  than  0.1  (Incropera  and  Dewitt, 
1985:179)  .  Consequently,  the  calculated  trace  for  the  second 
spacial  cell  should  match  an  exponential  temperature  rise  with 
a  thermal  time  constant  defined  by  T  =  p  where 
p  is  the  density,  is  the  thermal  capacity  of  the  steel,  V 
is  the  volume  of  steel,  h^onv  is  the  heat  transfer  coefficient, 
and  is  the  surface  area  of  the  matrix  (Ibid:  17 6)  .  A  test 


52 


case  was  performed  for  a  time  increment  At  =  10“^  seconds,  and 
twenty  one  space  increments.  The  comparison  is  shown  in 
Fig.  12.  The  calculated  trace  for  Tg^  matched  the  theoretical 
trace  with  a  mean  difference  of  3.1  x  10"®°  C  over  the  0.2 
seconds  of  the  transient.  Again,  the  comparison  is  very  good. 
Hence,  these  two  validation  tests  give  confidence  the 
numerical  model  is  working  as  intended. 

Determining  the  Heat  Transfer  Coefficient.  With  tsie 
available,  the  outlet  temperature  trace  can  be  obtained  for 
any  combination  of  geometric  and  flow  conditions,  and  a  value 
for  the  heat  transfer  coefficient.  The  actual  heat  transfer 
coefficient  can  be  determined  by  requiring  that  a 
characteristic  of  the  calculated  outlet  temperature  trace 
match  the  same  characteristic  of  the  experimental  temperature 
trace . 

Zeroin  Subroutine.  What  is  needed  then  is  a  way  to 
find  the  best  heat  transfer  coefficient.  This  particular 
problem  has  attracted  attention  for  centuries;  it  is  called  a 
root-solving  problem.  In  a  root-solving  problem,  one  seeks 
the  value  of  the  independent  variable  for  which  a  function  of 
the  variable  equals  zero,  in  the  current  research  fn(hcoNv)-0- 
Fortunately,  all  this  attention  has  resulted  in  an  algorithm 
by  Dekker  (1969),  improved  by  Brent  (1973)  as  reported  in 
Forsythe  (1977:  Section  7.2).  The  authors  claim  that  this 
algorithm,  called  Zeroin,  should  always  converge,  and  needs 
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less  than  eleven  iterations  to  converge  for  smooth  functions 
(Brent,  1973).  A  requirement  of  Zeroin  is  that  some  initial 
bounding  guesses  for  are  input  to  the  algorithm  such  that 

and  have  different  signs,  i.e.  one  is  greater 

than  zero  and  the  other  is  less  than  zero.  The  algorithm 
starts  by  evaluating  fnlh^.^^^)  for  the  initial  bounds  and 
corrects  this  value  until  the  best  heat  transfer  coefficient 
is  chosen.  This  is  determined  when  the  difference  between 
h^onv  f^rom  one  iteration  to  the  next  is  smaller  than  a  chosen 
tolerance.  For  the  current  research,  this  tolerance  was 
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chosen  to  be  <  10  0  W/m^/°C,  which  is  two  per  cent  of  the 
average  value  for  the  heat  transfer  coefficient  measured  in 
the  thirty  test  runs .  This  error  matches  the  error  in  the 
numerical  model  (1.8  %)  given  above  for  NMAX  =  100. 

The  main  Program.  A  FORTRAN  program  called  main 
controls  the  calculation  of  the  heat  transfer  coefficient  by 
employing  tsie  and  Zeroin  as  subroutines.  After  the  inputs  to 
main  are  made,  Zeroin  is  called,  which  in  turn  calls  tsie  to 
calculate  the  best  heat  transfer  coefficient  is  chosen 
as  described  above.  Main  also  calculates  the  friction  factor, 
the  effectiveness  of  the  regenerator,  and  the  theoretical  heat 
transfer  coefficient  based  on  the  approach  of  Pucci  et  al . 
(1967)  described  in  the  Theory  section.  These  other  results 
are  discussed  below. 

The  Three  Criteria.  Many  criteria  could  be  chosen 
from  which  the  heat  transfer  coefficient  is  determined.  The 
maximum  slope  of  the  output  temperature  trace  and  the  minimum 
root-mean- squared  difference  between  the  experimentally 
measured  rising  temperature  trace  and  an  analytical  solution 
were  mentioned  above.  In  the  current  research,  these  two 
criteria,  maximum  slope  and  minimum  root-mean-squared 
difference,  are  considered,  as  well  as  requiring  the  sponge 
effect  delay  time,  SEDT,  to  be  the  same  for  both  experimental 
and  calculated  results.  The  maximum  slope  and  delay  time 
cases  are  straightforward.  But  finding  the  heat  transfer 
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coefficient  which  minimizes  the  root-mean- squared  difference 
between  calculated  and  experimental  temperature  traces 
required  some  ingenuity. 

There  are  two  difficulties  defining  a  criterion  for  the 
minimum  root-mean-squared  difference  approach.  First,  a  root- 
mean-squared  difference  is  always  greater  than  zero.  As 
mentioned  above,  the  Zeroin  algorithm  requires  two  initial 
guesses  for  which  bound  its  value  and  for  which 

has  a  different  sign.  In  the  root-mean-squared  difference 
case,  fn(hco„„)  is  always  greater  than  zero.  To  get  around 
this  difficulty,  a  related  property  of  the  root-mean-squared 
difference  is  used.  Fig.  13  shows  the  anticipated  appearance 
of  the  root-mean- squared  difference  versus  h^^nv  curve.  The 
minimum  value  of  the  root-mean-squared  difference  is  where  the 
slope  is  zero.  The  slope  before  and  after  the  minimum  is  of 
different  signs.  Since  the  function  should  be  smooth,  the 
value  of  the  slope  at  any  h^onv  can  be  approximated  by  a  simple 
forward  difference,  i.e.  the  value  of  the  root-mean-squared 
difference  is  calculated  at  a  particular  value  of  h^onv 
also  at  a  value  some  small  amount  larger,  say  h^^nv  +  100.  The 
difference  in  these  two,  divided  by  the  difference  in 
100),  gives  the  local  value  of  the  slope.  Using  this 
approach,  Zeroin  and  tsie  can  be  used  to  solve  fnlh^o^^)  =  0 
where  the  function  is  the  slope  of  the  root-mean-squared 
difference  versus  h^onv  line.  When  the  slope  is  near  zero. 
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Shape  of  the  Heat  Transfer  Coefficient  vs  RMSD  Curve 


Figure  13.  Root -Mean-Squared-Difference  as  a  Function 

of  the  Heat  Transfer  Coefficient 


the  value  of  the  heat  transfer  coefficient  gives  the  best  fit 
to  the  numerical  data. 

The  second  problem  with  the  minimum  root-mean- squared 
difference  approach  is  deciding  over  which  time  interval  to 
compute  the  root-mean-squared  difference.  During  the  sponge 
effect  delay  time,  there  is  not  much  temperature  change  in  T,2- 
At  late  times,  the  physics  of  the  model  no  longer  fit  the 
physics  of  the  experiment  as  explained  in  the  Theory  section. 
In  order  to  avoid  some  amount  of  arbitrariness,  the  root-mean- 
squared  difference  is  evaluated  from  the  end  of  the  delay  time 
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to  one  matrix  thermal  time  constant  or  about  50  milliseconds 


after  the  occurrence  of  the  maximum  slope.  This  concentrates 
the  comparison  to  a  time  interval  which  includes  the  maximum 
slope,  i.e.  most  of  the  change  in  has  taken  place.  This 
completes  the  description  of  the  three  criteria.  A  brief 
description  of  other  results  calculated  by  main  concludes  this 
chapter . 


Other  Results.  Three  other  important  results  are 
returned  from  the  main  program:  1)  the  friction  factor,  2) 
the  effectiveness  of  the  regenerator,  and  3)  the  theoretical 
heat  transfer  based  on  the  maximum  slope  technique  of  Pucci, 
et  al .  (1967)  described  above. 

The  friction  factor  can  be  defined  in  many  ways  (Kays  and 
London,  1984:36;  Miyabe,  et  al ., 1982 : 1839 )  .  The  one  chosen 
for  this  research  is  developed  by  Armour  and  Cannon  (1968) 
specially  for  woven  screens . 


where  AP 

O 

PS 

L 

P 

and  Ug 


f 


AP  PS 
L  P 


(23) 


is  the  pressure  drop  in  the 
regenerator  (Pa) 
is  the  porosity 
is  the  matrix  pore  size  (m) 
is  a  length  =  Q  EL  where  Q  is  a  factor 
which  compensates  for  the 
"  tortuousity"  and  EL  is  the  matrix 
length  (m) 

is  the  gas  density  (kg/m^) 
is  the  gas  velocity  (m/sec) 
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The  AP  was  measured  for  each  test  run.  The  porosity  was 


determined  by  the  following  equation: 


=  1.0 


^su 
Pm  ^ 


(24) 


where 

^sa 

is 

the 

mass  of  the  SU  (kg) 

Pm 

is 

the 

density  of  the  matrix 

(kg/m') 

A 

is 

the 

tube  cross-sectional 

area  (m' 

and 

EL 

is 

the 

length  of  the  matrix 

(m) 

The  pore  size  was  determined  to  be 

PS={{{1/MESH) {2.0  RF))^/^  (25) 

where  MESH  is  the  mesh  size  (m’^) 

d„  is  the  wire  diameter  (m) 
and  RF  is  the  reduction  factor 

According  to  Armour  and  Cannon  (1968),  for  plain  weave 

screens,  the  length  should  be  the  length  of  the  matrix,  EL, 

and  Q  is  set  equal  to  1.0.  The  gas  density  is  calculated  from 

a  measured  average  temperature  and  pressure  in  the  test 

section,  and  perfect  gas  relations.  Finally,  the  mass 

velocity  is  the  measured  mass  flow  rate  divided  by  the  area 

for  flow,  Af  =  O  A,  divided  by  the  density. 

Another  parameter,  effectiveness,  was  defined  for  the 

purpose  of  deciding  the  merits  of  reducing  the  thickness  of 

the  screens.  Effectiveness  is  defined  for  regenerators  in 

Chapter  II  in  terms  of  the  ability  of  the  regenerator  to 
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remove  the  maximum  amount  of  energy  from  the  gas  relative  to 
an  ideal  case.  Another  effectiveness,  SPEFF,  can  be  defined 
for  the  case  of  one-directional  flow: 


SPEFF  = 


m 


Cp  SEDT 
^sn 


(26) 


n 


where  m  is  the  mass  flow  rate  (kg/ sec) 

Cp  is  the  specific  heat  of  the  gas  (J/kg/°C) 
SEDT  is  the  sponge  effect  delay  time  (sec) 

W^u  is  the  mass  of  matrix  (kg) 
and  Ci^  is  the  specific  heat  of  the  matrix 

(J/kg/°C) 

The  numerator  is  the  energy  of  the  gas,  with  respect  to  the 
initial  temperature  of  the  regenerator,  which  flows  into  the 
matrix  during  the  delay  time.  The  denominator  is  the  amount 
of  energy  required  to  raise  the  matrix  from  its  initial 
temperature  to  the  same  final  temperature  defined  by  the  flow 
of  the  gas.  If  the  regenerator  uses  up  all  of  its  thermal 
capacity  to  absorb  all  of  the  relative  energy  from  the  gas 
during  the  delay  time,  the  effectiveness  is  1.0.  More  will  be 
said  about  these  three  results  in  Chapter  V. 
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IV.  Experimental  Apparatus  and  Procedure 


This  chapter  describes  the  apparatus  and  procedures  used 
to  gather  the  experimental  data,  and  the  data  reduction 
technique  for  composing  the  input  data  file  for  each  test  run. 

Apparatus 

First,  details  of  the  experimental  apparatus  are  given. 
A  diagram  is  shown  in  Fig.  14.  A  detailed  list  of  the 
equipment  for  Fig.  14  is  given  in  Table  1. 


Figure  14 .  Diagram  of  Eaq>erimeiital  Apparatus 
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Table  1.  Instrumentation  List  for  Fig.  10 


ITEM  # 

DESCRIPTION 

SERIAL  # 

1 

Chiller  -  NESLAB  HX-lOO 

A95064014 

2 

4  Tank  Helium  Harness  -  Cold  Gas 

-- 

3 

Single  Helium  Tank  -  Hot  Gas 

-- 

4 

Regulator  -  Victor  SR  -  600 

DB  48406 

5 

Same  as  Item  #4 

DB  48407 

6 

Heating  Coil  -  Hot  Water  Bath 

-- 

7 

Heater  -  Fisher  Thermix  11-493 

121  MX  2614 

8 

Heat  Exchanger  -  Counter  Flow 

1 

9 

4  Way  Valve  -  Whitey  B-45YF8 

-- 

10 

2  X  Thermistor  -  Thermometrix  FP07 

2  X  Motherboard  -  E  &  L  Instruments 

-- 

11 

Regenerator 

-- 

12 

Flowmeter  -  Venturi  Tube 

-- 

13 

4  X  Valve  -  Whitey  B-44S6 

-- 

14 

Pressure  Transducer  -  Validyne  DPlO-42 

74693 

15 

Pressure  Transducer  -  Validyne  DP15-48 
"  -  ”  DP15-56 

87548 

84654 

16 

Data  Acquisition  System  -  Zenith  Z-510 

4LSBUX7536 

17 

Thermocouple  -  Omega  .005"  Fe-constan 
Indicator  -  Omega  DP41-TC-A 

10984A1-01 

The  first  major  subsection  of  the  apparatus  is  the 
equipment  used  to  handle  the  helium.  Helium  was  chosen  as  the 
working  fluid  since  99.995%  pure  helium  was  available  (small 
pores  in  a  matrix  are  easily  clogged) ,  and  helium  is  commonly 
used  for  the  working  fluid  in  cryocoolers,  although  any  clean 
gas  could  have  been  used.  The  tanks  used  to  hold  the  helium 


62 


gas  have  17  MPa,  0 . 5  capacity  and  hold  about  5  kg.  A  four- 
tank  harness  supplied  the  cold  gas  flow.  A  single  tank  was 
used  for  the  hot  gas  flow.  A  NESLAB  Model  HX-lOO  chiller 
cooled  water  to  3°  C.  The  water  was  pumped  into  a  3.0  m 
length  section  of  25  cm  diameter  PVC  piping.  This  forms  the 
outer  shell  of  a  counter- flow  heat  exchanger.  The  helium 
flowed  through  approximately  50  meters  of  coiled  copper  tubing 
over  which  cold  water  was  pumped.  The  gas  could  be  cooled  to 
6°  C.  A  Fisher  Thermix  adjustable-power  hot  plate  was  used  to 
heat  a  steel  container  filled  with  water.  The  hot  gas  flowed 
through  a  3  m  length  of  9.5  mm  O.D.  tubing,  coiled  25  cm  in 
diameter,  which  lay  in  the  hot  water.  The  gas  could  be  heated 
to  50°  C.  The  gas  flowed  through  the  test  section  due  to  the 
pressure  in  the  tanks  which  was  regulated  by  Victor  SR- 600 
regulators.  The  mass  flow  rate  was  controlled  with  tandem 
upstream  and  downstream  needle  valves.  Except  for  the  tubing 
of  the  SU  in  the  test  section,  the  gas  flowed  through  9.5  mm 
O.D.  copper  tubing  rated  to  34.5  MPa.  The  tubing,  valves, 
geometry,  etc.,  were  chosen  to  achieve  a  Reynolds  number 
(based  on  pore  size)  between  50  and  600. 

The  next  major  subsection  of  the  experimental  apparatus 
is  the  test  section.  At  the  test  section  entrance  is  a  Whitey 
Model  B-45YF8  four-way  valve  (Fig.  14) .  With  the  valve  handle 
in  the  vertical  position,  cold  gas  flowed  through  the  test 
section;  in  the  horizontal  position,  hot  gas  flowed  through 
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the  test  section.  The  flow  diverted  from  the  test  section  was 


vented  to  atmosphere;  the  outflow  from  the  test  section  was 
also  vented  to  atmosphere.  The  central  element  of  the  test 
section  is  the  fixture  for  holding  the  screen  units  (SUs)  and 
instrumentation  for  pressure  and  temperature  measurements. 

One  of  the  two  mounting  fixtures  on  either  end  of  an  SU 
is  shown  in  Fig.  15.  They  were  made  from  0.0095  m  I.D.  to 
0.0159  m  I.D.  Swagelok  brand  male  connectors.  The  connectors 
were  reamed  such  that  the  SUs  slid  far  enough  down  the 
interior  of  the  opening  that  the  Swagelok  gas-tight  seal  would 
clamp  approximately  0.0032  m  from  the  end  of  the  SU.  Three  100 
mesh  screens  were  inserted  into  the  upstream  connector  inlet 
to  act  as  a  mixer  and  turbulence-generator.  Each  connector 
held  three  measurement  devices:  1)  a  Thermometrix  Fas tip  1.5 
X  10'^  m  diameter  bead,  quick-response,  thermistor,  2)  a 
0.00635  m  O.D.  f lushed-wall ,  copper  tube  for  pressure 
measurements,  and  3)  a  TSI  brand  hot-wire  probe  (not  used  for 
measurements  in  this  report) .  The  three  measurement 
instruments  were  epoxied  into  the  connectors  with  specially 
made  holding  fixtures.  The  connectors  were  slid  over  the  end 
of  the  SU  to  a  predetermined  position,  and  then  clamped  in  a 
leak-proof  seal  with  Swagelok  fittings.  The  SUs  were  easily 
removed  and  replaced  for  testing. 

The  design  of  the  SUs  is  shown  in  Fig.  16.  The  SUs  are 
made  from  304  stainless  steel  plain-weave  mesh  in  three  sizes: 
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Figure  15.  Instrumentation  Holders 


100,180,  and  250  strands  per  inch.  The  screen  was  rolled  by 
a  collandering  process  (between  a  stationary  rotating  pin  and 
a  moving  table)  on  a  Bridgeport  mill  to  nominal  15%,  30%  and 
50%  reduction  in  thickness,  i.e.  reduction  factors  of  0.85, 
0.70,  and  0.50.  Unrolled  screens  were  also  tested.  A  special 
tool  was  used  to  punch  out  disks  of  a  diameter  that  fit  snugly 
inside  the  stainless  steel  tube  used  to  hold  the  regenerator 
matrix.  There  are  precisely  two  hundred  screens  per  SU  since 


Figure  16.  Photo  of  Screen  Unit  (S/U) 


this  allows  direct  comparison  of  SU  performance.  The  number 
two  hundred  was  chosen  because  studies  show  (Mikulin,  et  al .  , 
1972)  that  this  is  a  sufficient  number  of  screens  to  eliminate 
length  effects  on  heat  transfer.  The  screens  were  cleaned 
with  isopropyl  alcohol,  dried,  and  stacked  inside  a  0.0159  m 
O.D.,  5.1  X  10“^  m  thick  walled,  ASTM-A269  304  stainless  steel 
tube.  The  weights  and  lengths  of  the  individual  screens  and 
entire  SUs  were  carefully  measured  with  vernier  calipers  and 
a  Denver  Instruments  XL-300  electronic  scale.  A  9.5  mm  long 
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spacer  made  from  PVC  tubing  was  epoxied  into  each  end  of  the 
tube  to  hold  the  screens  in  place. 

Thirteen  SUs  were  fabricated  as  listed  in  Table  2 .  The 
first  twelve  are  from  the  three  mesh 'sizes  and  four  reduction 
factors. 


Table  2 .  List  of  Screen  Units 


s/ 

MESH 

Porosity 

RF^- 

dj 

p^c. 

EL 

d. 

ts 

e. 

^su- 

u 

3) 

10- 

®m 

10-^m 

10- 

'm 

10' 

®m 

gms 

n 

100 

0 

.  677 

1.0 

11 

.4 

16.5 

45 

.4 

22 

.  6 

19.0 

2 

M 

0 

.  636 

0.86 

II 

15 . 6 

38 

.9 

19 

.  6 

19.1 

3 

II 

0 

.599 

0.70 

II 

14.6 

31 

.9 

15 

.9 

18.8 

4 

II 

0 

.466 

0.50 

II 

13.1 

22 

.7 

11 

.2 

17.5 

5 

180 

0 

.660 

1.0 

6. 

4 

9.15 

29 

.5 

14 

.7 

12.1 

6 

11 

0 

.601 

0.83 

II 

8.60 

24 

.5 

12 

.3 

12.1 

7 

II 

0 

.520 

0.72 

II 

8.20 

21 

.2 

10 

.5 

12.0 

8 

II 

0 

.375 

0.53 

II 

7.41 

15 

.7 

7  . 

4 

11.6 

9 

250 

0 

.  666 

1.0 

9 

6.71 

17 

.4 

8. 

6 

6.9 

10 

II 

0 

.598 

0.85 

II 

6.35 

14 

.7 

7  . 

4 

6.6 

11 

II 

0 

.541 

0.75 

M 

6.08 

13 

0 

6. 

5 

6.4 

12 

II 

0 

.399 

0.55 

II 

5.51 

9. 

6 

4. 

7 

6.3 

13 

MIX 

0 

.577 

MIX 

7. 

3 

9.28 

23 

.1 

11 

.8 

12.5 

a.  Reduction  Factor;  b.  Wire  diameter;  c.  Pore  Size;  d.  Length;  e. 
Screen  thickness;  f.  Mass  of  the  SU 


A  special  SU  was  also  made  to  investigate  an  alternate 
configuration  which  may  achieve  the  desired  result,  i.e.  ,  less 
void  fraction  with  comparable  compactness  factor.  The  special 
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SU  was  made  from  a  random  mixture  of  all  the  mesh  sizes  and 


reduction  factors . 

Pressure  measurements  were  made  at  either  end  of  the  SU 
to  determine  the  pressure  drop,  and  subsequently,  a  friction 
factor.  The  pressure  at  the  inlet.  Pi,  was  measured  with  a 
Validyne  Model  DP-1556  34.5  MPa  differential  pressure 
transducer.  The  differential  pressure  transducer  measured  the 
difference  between  PI  and  atmospheric  pressure  which  was 
monitored  by  a  Transamerica  Delaval  Type  CEC  2500  Digital 
Barometer.  The  pressure  drop  across  the  SU  was  measured  with 
a  Validyne  Model  DP-1548  0.5  MPa  differential  transducer. 
Both  transducers  used  a  Validyne  Model  CD-280  six  channel 
exciter . 

Temperature  measurements  were  made  in  three  locations . 
At  the  inlet  and  outlet  from  the  SU,  a  thermistor  was  used. 
The  thermistor  was  one  leg  in  a  Wheatstone  bridge 
configuration.  A  Hewlett-Packard  Model  6405C  DC  Power  Supply 
provided  1.0  ±  0.001  VDC  across  the  bridge.  This  low  power 
voltage  was  used  to  minimize  self-heating  of  the  thermistors. 
The  major  benefit  of  the  bridge  was  that  the  other  three 
resistances  could  be  chosen  so  that  the  bridge  operated  near 
a  balanced  condition.  The  smaller  the  voltage  measurement 
range  for  the  thermistor  (±  0.1  VDC  was  achieved  in  this 
case) ,  the  greater  the  amplification  by  the  data  acquisition 
system  (DAS) ,  and  the  better  the  digital  resolution  of  the 
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fluctuating  component  of  the  signal  (Labview  for  Windows, 
1993;  2-8) .  The  third  location  for  temperature  measurements 
was  at  the  entrance  to  the  bypass  vent.  The  temperature  there 
was  measured  by  a  0.127  mm  (0.005")  diameter  Omega  iron- 
constantan  thermocouple  in  conjunction  with  an  Omega  Model 
DP41-TC-A  Digital  Thermocouple  Indicator.  The  measurement  was 
made  to  monitor  the  bypass  gas  temperature,  in  particular,  to 
determine  when  it  was  hot  enough  to  begin  the  transient  of  a 
test  run. 

At  the  exit  of  the  test  section,  a  venturi  flowmeter  with 
a  0.0065  m  diameter  throat  was  used  for  measuring  the  mass 
flow  rate  of  the  gas.  The  venturi  flowmeter  had  geometrical 
proportions  recommended  by  ASME  standards  (ASME,  1971) .  The 
pressure  drop  across  the  throat  was  measured  with  a  Validyne 
Model  DP  10-42  34.5  kPa  differential  pressure  transducer  using 
the  same  exciter  mentioned  above. 

The  last  major  subsection  of  the  experimental  set-up  is 
the  Data  Acquisition  System  (DAS) .  Coaxial  cables  were  used 
to  connect  the  instrument  readings  (all  DC  voltages)  to  a 
junction  board.  The  board  was  connected  to  a  Zenith  Data 
Systems  Z-Station  510  with  a  ribbon  cable.  All  amplification 
and  conversions  were  done  internally  to  the  computer,  using  a 
National  Instruments  NI-DAQ  driver  and  an  AT-MIO-16  board 
which  is  controlled  by  Labview  Version  3.0.1  Virtual 
Instrument  software.  Special  Labview  Virtual  Instruments  (VI) 
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were  designed  and  built  for  calibrations  and  experimental 
control.  MY_SS_REGENERATOR_HTC.VI  provided  steady-state 
monitoring  of  conditions  in  the  test  section  and  transient 
data  collection,  conversion,  and  storage.  Once  the  data  were 
stored  to  a  file,  they  could  be  transferred  via  the  local  area 
network  to  a  SUN  workstation  where  data  reduction  was 
accomplished. 

Experimental  Procedure 

The  detailed  experimental  procedure  is  shown  in  Appendix 
A.  A  synopsis  is  given  here.  The  preliminary  part  of  the 
procedure  involved  preparing  for  the  test.  After  a  new  SU  was 
placed  into  the  test  section,  a  leak  test  was  done  around  all 
seals .  Foam  insulation  was  placed  on  any  surface  open  to  the 
room  air.  Voltages  were  checked  on  the  thermistor  circuits 
with  a  Hewlett-Packard  Model  34401A  Multimeter.  Gas  levels 
were  maintained  high  enough  to  complete  the  test,  and  the 
chiller  and  heating  units  were  brought  to  their  respective 
steady  state  conditions. 

The  testing  procedure  began  with  a  data  sheet  containing 
the  appropriate  information  for  the  test  run  about  to  begin. 
The  approximate  back  pressure  was  set  in  the  tanks  and  both 
hot  and  cold  flows  were  adjusted  so  that  the  correct  pressure. 
Pi,  at  the  entrance  to  the  SU,  and  the  correct  mass  flow  rate 
were  maintained.  After  approximately  four  minutes,  steady 
state  cold  temperatures  around  10°  C  prevailed  in  the  SU,  and 
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hot  gas  temperatures  around  35°  C  prevailed  in  the  vent. 
Next,  the  data  collection  trigger  was  tripped,  and  after  an 
approximately  three  second  delay,  the  four-way  valve  was 
manually  rotated  so  that  the  hot  gas  was  diverted  to  the  SU. 
Pre- transient  conditions  were  recorded  by  the  DAS,  and  data 
were  taken  during  the  transient.  After  approximately  20 
seconds,  all  gas  flows  were  shut  off.  A  graph  of  the  test  run 
data  was  shown  on  the  computer  screen  so  that  the  quality  of 
the  test  run  could  be  ascertained.  Post-transient  data  as 
well  as  any  pertinent  comments  were  added  to  the  data  sheet 
for  the  test.  The  data  file  was  given  an  identifier,  backed- 
up  on  tape,  and  sent  to  a  SUN  workstation  via  a  file  transfer 
protocol  (ftp)  and  the  local  network.  The  data  reduction 
technique  described  in  the  next  section  was  applied  to  the 
data  before  it  was  loaded  into  the  numerical  model  for 
determination  of  heat  transfer  and  friction  coefficients. 

Data  Reduction  Technique 

After  temperature,  pressure,  and  mass  flow  data  had  been 
gathered  for  the  transient  tests,  they  were  converted  into  a 
format  suitable  to  determine  the  heat  transfer  coefficient. 
The  aim  of  this  section  is  to  describe  the  steps  taken  to 
create  a  data  file  that  was  fed  into  the  FORTRAN  program  which 
models  the  flow.  Some  important  parameters  are  defined,  as 
well  as  how  they  were  calculated  from  the  experimental  data. 
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The  raw  data  from  the  data  acquisition  system  was  sent  to 
a  Sun  SparclO  workstations  via  file  transfer  protocol  over  the 
local  area  network.  This  data  contained  temperature  traces 
from  the  inlet  and  outlet  to  the  test  section,  pressure  and 
pressure  drop  readings  from  the  three  pressure  stations  in  the 
test  section,  and  an  average  value  of  the  mass  flow  rate.  The 
data  sampling  rate  and  total  transient  time  were  also 
included.  These  data  were  loaded  into  a  MATLAB  session  for 
manipulation.  MATLAB  is  a  high-performance  numerical 
computation  and  visualization  software  (MATLAB,  1993).  The 
MATLAB  m. files  used  to  perform  the  data  reduction  are 
contained  in  Appendix  B . 

The  first  important  parameters  determined  from  the  data 
were  the  start  and  finish  times  of  the  transient.  Each  test 
run  lasted  ten  seconds  and  data  were  collected  at  10,000  time 
steps  (a  sampling  frequency  of  1  kHz) .  To  is  the  start  time 
and  Xf  the  finish  time  of  the  transient  which  began 
approximately  three  seconds  after  data  acquisition  began. 
These  data  were  used  to  concentrate  the  matching  of  the 
analytical  and  experimental  results  during  the  transient.  The 
raw  data  for  and  were  filtered  with  a  20th  order  finite 
impulse  response  (FIR)  filter.  Xq  was  defined  to  be  the  time 
at  which  changed  by  more  than  two  standard  deviations  (0.4° 
C)  from  its  running  average  prior  to  the  transient.  Xf  was 
chosen  more  subjectively.  The  MATLAB  m.file  (trange.m). 
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written  to  choose  Tf  starts  by  identifying  the  data 
acquisition  time  that  the  slope  of  the  trace  is  greater 
than  five  degrees  per  second.  This  normally  gave  a  long 
duration  value  for  the  transient  duration.  The  value  of  Xf 
was  shortened  from  the  MATLAB  calculation  to  a  point  after 
which  the  maximum  slope  in  the  outlet  temperature  trace  had 
obviously  occurred,  between  ten  and  fifteen  matrix  time 
constants  after  Xq. 

Once  the  beginning  of  the  change  in  the  inlet  temperature 
trace  and  end  of  the  change  in  the  outlet  temperature  trace 
were  identified,  the  next  important  parameter  to  find  was  Xj . 
The  time  between  Xq  and  when  the  trace  of  begins  to  change 
measurably,  X2,  is  defined  as  the  sponge  effect  delay  time, 
SEDT.  The  determination  of  X2  was  done  in  the  same  way  as  for 
Xq,-  by  definition,  it  s  located  at  the  data  acquisition  time 
step  where  Tg2  has  changed  by  more  than  two  standard  deviations 
from  its  running  average  since  the  test  run  began. 

The  final  important  parameters  were  the  magnitude  of  the 
maximum  slope  of  the  T^2  trace  and  the  time  step  at  which  it 
occurred.  A  MATLAB  m.file  (bestpoly.m)  was  constructed  to 
find  the  best  polynomial  fit  to  the  trace  during  the  time 
between  X2  and  Xf.  It  picked  the  polynomial  of  order  two 
through  thirty  which  best  fits  the  experimental  data  in  a 
root-mean-squared-difference  sense.  Another  m.file 
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(maxslope .m)  calculated  the  maximum  slope  of  the  polynomial, 
and  when  it  occurred  in  relation  to  Xq. 

The  final  step  in  the  reduction  process  was  the 
collection  of  all  the  data  into  a  read  file  for  each 
experimental  run.  Results  from  the  raw  data  such  as  mass  flow 
rate  and  pressures  were  assembled  at  the  beginning  of  the  read 
file.  Next,  some  geometric  properties  of  the  individual  SUs 
such  as  mass,  pore  size,  and  length  were  added  to  the  file. 
Next,  the  filtered  and  traces  were  appended  to  the 
file.  Finally,  the  x  parameters  described  above  were 
included.  All  these  data  were  used  by  the  numerical  model 
given  in  the  next  section  to  determine  a  heat  transfer 
coefficient . 
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V.  Results  and  Discussion 


In  this  chapter,  results  which  give  confidence  in  the 
method  described  in  Chapter  III,  and  which  prove  the 
hypothesis  stated  in  Chapter  I,  are  presented.  Additionally, 
a  section  is  dedicated  to  comparisons  between  the  heat 
transfer  coefficients  obtained  using  the  three  criteria 
mentioned  in  Chapter  III.  How  regenerator  effectiveness  is 
changed  by  reducing  the  screen  thickness  is  shown,  and  a 
general  discussion  about  the  results  is  given  at  the  end  of 
the  chapter.  The  numerical  values  for  the  test  results  are 
listed  in  Appendix  E.  Before  the  test  results  are  shown,  some 
preliminary  results  and  definitions  are  presented. 


Preliminary  Results 

Four  preliminary  items  are  addressed.  One  is  a 
definition  of  the  Reynolds  number;  the  second  is  a  relation 
for  specific  surface  area  for  each  matrix;  the  third  is  a 
typical  temperature  trace  from  the  experimental  apparatus;  and 
the  fourth  is  the  test  run  matrix. 

During  the  course  of  this  research,  many  definitions  for 
the  Reynolds  number  were  encountered.  Generally,  it  is 
defined 


Re 


P  UgL 

P 


(27) 
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where 


Re  is  the  Reynolds  number 

p  is  the  gas  density  (kg/m^) 

Ug  is  the  gas  velocity  (m/sec) 

|i  is  the  dynamic  viscosity  (kg/m/ sec) 
and  L  is  a  characteristic  length  (m) 

The  difference  between  Reynolds  numbers  is  the  choice  of  the 
characteristic  length.  The  most  common  characteristic  length 
in  the  porous  media  arena  is  the  pore  size.  The  pore  size, 
PS,  is  defined  as  the  cube  root  of  the  average  volume  of  the 
voids  in  the  medium.  In  the  study  of  heat  transfer  and  fluid 
dynamics  in  pipes  and  over  cylinders,  it  is  also  common  to  see 
a  Reynolds  number  based  on  wire  diameter  or  hydraulic  radius . 
For  a  stacked-screen  matrix,  the  choice  is  not  clear.  In  an 
article  written  by  Armour  and  Cannon  (1968) ,  a  definition  is 
developed  specially  for  the  case  of  stacked  screen  matrices. 
They  contend  the  pressure  drop  in  a  matrix  is  the  sum  of  two 
terms:  1)  a  surface  shear /drag  term  which  requires  an  area, 
and  2 )  an  inertial  term  due  to  eddies  and  sudden 
enlargements/contractions  in  the  direction  of  flow  which 
requires  a  volumetric  dimension.  The  resulting  definition  for 
the  characteristic  length  is 


L 


1 . 0 


'A  p 

aI  PS 


(28) 
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where 


is  the  surface  area  per  unit  volume  of 
unrolled  matrix  (m”^) 
and  PS  is  the  pore  size  (m) 

This  combination  of  area  and  linear  dimensions  was  shown  by 

Armour  and  Cannon  to  give  the  best  correlations  for  conditions 

in  a  stacked-screen  matrix  and  will  be  used  in  the  current 

research  unless  otherwise  stated. 

In  order  to  use  the  Aormour  and  Cannon  dimension,  the 

surface  area  per  unit  volume  needs  to  be  determined  for  each 

matrix.  This  was  accomplished  using  geometrical  arguments  and 

measurements  of  the  dimple  in  the  wire  left  by  the  rolling 

process.  Wiese  (1993:4-2)  took  electron  microscope 

photographs  of  the  screens  which  he  rolled  by  the  same  process 

(same  machine)  to  the  same  reduction  factors  as  those  used  in 

the  current  research.  The  photographs  show  that  the  rolling 

caused  the  round  wire  surfaces  to  be  flattened  at  certain 

locations  in  each  screen  cell  (set  of  four  crossing  wires) . 

When  two  screen  disks  are  stacked  against  one  another,  the 

flattened  areas  cover  each  other,  and  these  covered  areas  no 

longer  come  in  contact  with  the  fluid.  Wiese's  photographs 

were  used  to  estimate  the  area  per  unit  volume  for  each  of  the 

rolled  screen  matrices.  The  results  and  a  simple  curve-fit 

are  shown  in  Fig.  17. 

To  use  this  plot,  the  specific  area  per  unit  volume  for 
unrolled  screens,  A^,,  is  obtained  from  a  geometrical  relation 
developed  in  Armour  and  Cannon  (1968:418). 
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where  PIT  is  the  inverse  of  the  mesh  number  (m) 

and  d„  is  the  wire  diameter  (m) 

This  value  for  A^,  is  multiplied  by  the  appropriate  factor 

shown  in  Fig.  17  to  get  the  specific  area  of  the  rolled  screen 

per  unit  of  unrolled  volume, 
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The  third  piece  of  preliminary  information  is  a  plot  of 
a  typical  inlet  and  outlet  temperature  trace.  This  is  shown 
in  Fig.  18.  The  temperatures  of  the  inlet  and  outlet  gas 
begin  at  steady  state  at  around  10°  C,  and  rise  quickly  to 
between  3  0°  C  and  3  5°  C.  The  plot  shows  the  filtered  data 
which  is  input  to  the  numerical  model. 

The  final  piece  of  preliminary  information  is  the  test 
run  matrix.  This  is  shown  in  Table  3.  Four  pieces  of 
information  are  given:  1)  the  test  run  number,  2)  the  SU 


Test  Run  #r10t034 


01  23456789  10 


_ Test  Run  Time  -  Sec _ 

Figure  18.  Typical  Eiiperimental  Temperature  Trace 
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number  (the  characteristics  for  which  can  be  found  in  Table 
2),  3)  the  Reynolds  number  based  on  pore  size  which  is  more 
often  used  in  the  literature  to  report  operating  conditions, 
and  4)  the  nominal  mass  flow  rate  for  each  test.  The 
geometries  and  Reynolds  numbers  are  representative  of  current 
operating  conditions  for  regenerative  refrigeration  cycles . 

Comparisons  to  Published  Data 

In  this  section,  published  data  for  unrolled  screens  of 
similar  geometry  and  flow  conditions  are  compared  to  the 
experimental  data  collected  during  the  current  research. 

In  Fig.  19,  the  data  for  the  three  SUs  made  from  unrolled 
screens  is  compared  to  published  correlations  for  friction 
factors  and  Nusselt  numbers.  In  the  case  of  the  friction 
factors.  Armour  and  Cannons's  (1968)  result,  /=  8.61/Re  +0.52, 
was  used  (solid  line) .  Their  study  included  similar  mesh 
sizes  and  a  Reynolds  number  range,  1.0  <  Re  <  100.  The 

agreement  is  very  good  with  a  mean  difference  of  ±  2.5%. 

The  published  Nusselt  number  data  (solid  line)  are  due  to 
Hamaguchi  (1990),  Nu  =  0.42  Re^^”'^®.  He  also  used  identical 
mesh  sizes,  and  a  Reynolds  number  range  of  1.0  <  Re  <  900. 
The  fit  is  also  good  with  a  mean  difference  of  -13.9%.  For 
the  reasons  given  in  Chapter  III,  i.e.  non-exact  step-change 
in  inlet  temperature  and  the  thermal  inertia  of  the  tube,  heat 
transfer  results  from  the  current  research  should  be  larger 
than  previously  published  ones.  These  favorable  comparisons 
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Table  3 .  Test  Run  Matrix 


Test  Run  # 

Screen  Unit  # 

Reynolds  number 
(Pore  Size) 

Mass  Flow  Rate 
(g/sec) 

01 

09 

50 

1.87 

02 

09 

100 

3.74 

03 

10 

50 

1.87 

04 

10 

100 

3.74 

05 

11 

50 

1.87 

05 

11 

50 

1.87 

06 

11 

100 

3.74 

07 

12 

50 

1.87 

08 

12 

100 

3.74 

09 

01 

250 

3.78 

10 

01 

500 

7.55 

11 

02 

250 

3.78 

12 

02 

500 

7.55 

13 

03 

250 

3.78 

13 

03 

250 

3.78 

14 

03 

500 

7.55 

15 

04 

250 

3.78 

16 

04 

500 

7.55 

17 

05 

100 

2.72 

17 

05 

100 

2.72 

18 

05 

250 

6.80 

19 

06 

100 

2.72 

20 

06 

250 

6.80 

21 

07 

100 

2.72 

22 

07 

250 

6.80 

23 

08 

100 

2.72 

23 

08 

100 

2.72 

24 

08 

250 

6.80 

25 

13 

150 

2.72 

26 

13 

200 

3.78 

Reynolds  Number  based  on  Wire  Diameter  -  Redw 


Figure  19.  Friction  and  Nusselt  Number  Comparisons 


give  confidence  that  the  method  outlined  in  Chapter  III  is 
acceptable  and  gives  believable  results  for  the  heat  transfer 
coefficient  and  friction  factors  of  the  unrolled  screen 
matrices,  and  will  also  work  for  the  rolled  screen  matrices. 

Repeatability 

Another  test  of  confidence  for  the  methodology  is 
repeatability.  In  Fig.  20,  data  is  shown  for  both  friction 
factor  and  Nusselt  number  for  four  sets  of  test  runs. 
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Demonstrated  Repeatability  of  Nusselt  Number  and  Friction  Factor 


Figure  20.  Repeatability  of  the  Data 

The  test  runs  were  done  on  different  days  but  with 
approximately  the  same  operating  conditions .  The  test  runs 
were  chosen  at  random  and  cover  a  range  of  reduction  factors 
and  Reynolds  numbers .  The  friction  data  is  particularly  good 
with  a  mean  difference  of  only  2.0%  between  the  members  of 
each  set.  The  fit  for  the  Nusselt  numbers  is  also  good  with 
a  mean  difference  of  8.4%  between  members  of  the  same  set. 
Since  the  results  can  be  duplicated  so  closely  from  one  test 
run  to  the  next,  confidence  in  the  measurements  was  warranted. 
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Heat  Transfer  Coefficient 


In  this  section,  the  results  for  the  heat  transfer 
coefficient  are  shown.  The  Nusselt  Number  is  defined  as 
Nu  =  where  h^onv  is  the  convective  heat  transfer 

coefficient,  is  the  Armour  and  Cannon  (1968)  dimension,  and 
k  is  the  thermal  conductivity  of  the  gas .  The  general 
relation  is  shown  in  Fig.  21. 
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The  best  fit  for  the  data  is  a  familiar  form  for  heat  transfer 
from  a  surface  to  a  fluid,  Nu  =  A  *  Re®,  where  in  this  case 
A=0.21  ±0.049  and  B=0.62  ±0.055. 

One  thing  seems  to  be  missing  from  the  relation  for  the 
Nusselt  number.  There  is  no  dependence  on  the  reduction 
factor.  During  the  first  attempt  to  run  the  data  in  the 
numerical  model,  the  surface  area  of  the  regenerator  was 
assumed  to  be  the  same  for  the  reduced  thickness  screens  as  it 
was  for  the  unrolled  screens.  When  this  was  the  case,  the 
Nusselt  number  decreased  as  the  reduction  factors  decreased. 
When  the  reduction  in  area  caused  by  rolling  the  screens  was 
taken  into  account,  by  using  a  characteristic  dimension  which 
included  the  effects  of  the  surface  area,  no  dependence  on  the 
reduction  factor  by  the  heat  transfer  coefficient  was 
observed.  It  appears  no  new  heat  transfer  mechanism,  e.g., 
increased  turbulence  intensity  or  additional  turbulence  scale, 
is  introduced  by  rolling  the  screens. 

In  Fig.  22,  the  reduced  Nusselt  number,  defined  as  the 
actual  Nusselt  number  divided  by  the  Nusselt  number  calculated 
for  an  unrolled  screen,  is  shown  plotted  against  the  reduction 
factor.  The  plot  shows  within  the  accuracy  of  the  testing, 
changing  the  reduction  factor  has  no  effect  on  the  heat 
transfer  except  by  changing  the  surface  area. 
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Reduced  Nusselt  Number  vs  Reduction  Factor 


Figure  22.  Reduced  Nusselt  Number  as  a  Function  of 

Reduction  Factor 


Friction  Factor 

In  this  section,  the  results  for  the  friction  factor  are 
given.  The  behavior  of  the  friction  factor  differs  from  the 
behavior  of  the  Nusselt  number.  The  reduction  factor  plays  a 
significant  role  in  the  value  of  the  pressure  drop. 

The  relation  between  the  friction  factor  and  Reynolds 
number  is  shown  in  Fig.  23.  The  plotted  line  is  the 
relationship  by  Armour  and  Cannon  (1968)  for  unrolled  screens. 
The  first  thing  to  notice  about  the  experimental  data  is  a 
group  of  points  which  are  far  from  the  rest  in  the  upper  left 
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hand  corner.  All  of  these  data  points  are  for  test  runs  with 
a  reduction  factor  =  0.5.  Apparently,  reducing  the  thickness 
of  the  screen  by  50%  causes  a  large  increase  in  the  friction 
coefficient . 

The  second  thing  to  notice  in  Fig.  23  is  the  remaining 
data  points  seem  to  fall  closely  to  the  line  for  the  unrolled 
relationship.  To  insure  this  was  the  case,  the  reduced 
friction  factor  was  plotted  against  the  reduction  factor  as 
shown  in  Fig.  24.  The  reduced  friction  factor  is  defined  as 
the  actual  friction  factor  divided  by  the  friction  factor  for 
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the  unrolled  screens  at  the  same  Reynolds  number.  It  appears 
that  the  friction  factor  increases  for  all  values  of  reduction 
factor  less  than  0.7. 

A  closer  look  at  the  friction  factor  results  for 
0.7  <  RF  <  1.0  is  shown  in  Fig.  25.  For  RF  =  0.7,  the 

friction  factor  increases  by  approximately  10%.  Although  the 
effect  of  increasing  the  friction  factor  is  small  in  this 
range,  it  is  nonetheless  present.  These  data  indicate  that 
the  friction  factor  increases  whenever  the  screens  are 
flattened  by  30%  or  more.  This  conclusion  is  consistent  with 
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the  results  shown  in  the  bottom  part  of  Fig  5.  More  will  be 
said  about  this  later. 

Compactness  Factor 

The  most  important  result  for  the  purpose  of  this 
research  is  the  compactness  factor.  The  compactness  factor  is 
a  ratio  between  the  heat  transfer  and  pressure  drop 
characteristics  of  a  regenerator.  The  definition  of 

compactness  factor:  CF  =  where  j„  is  the  Colburn  factor, 

and  /  is  the  friction  factor  defined  in  Chapter  III.  The 
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Colburn  factor  is  defined:  jn  =  St  where  St  is  the 

Stanton  number  and  Pr  is  the  Prandtl  number.  The  Stanton 
number  is  defined:  St  =  Nu  /(Re  Pr)  where  Nu  is  the  Nusselt 
number,  and  Re  is  the  Reynolds  number  defined  above. 

In  Fig.  26,  the  results  for  the  compactness  factor  for 
all  the  SUs  is  shown.  While  the  general  trend  is  lower  CF  for 
lower  values  of  the  reduction  factor,  there  is  no  discernable 
trend  which  depends  on  the  MESH  size. 

In  Fig.  27,  the  results  for  the  compactness  factor  as  a 
function  of  reduction  factor  are  shown  for  each  mesh  size. 
The  results  for  the  unrolled  cases  compare  favorably  with 
published  results  for  the  compactness  factor  (Radebaugh  and 
Louie,  1986:180).  For  the  100  and  180  mesh  screens,  the 
trend  is  clear,  the  compactness  factor  decreases  with 
decreasing  reduction  factor.  Hence,  rolling  the  screens  tends 
to  reduce  the  compactness  factor. 

For  the  250  mesh  screens,  the  compactness  factor  seems  to 
increase  at  a  reduction  factor  of  0.7.  The  data  for  this  case 
was  examined  more  closely,  and  whereas  the  heat  transfer 
coefficient  is  about  what  one  would  expect  (compared  to  the 
other  data) ,  the  friction  factor  is  smaller  than  expected. 
The  reason  is  not  apparent  for  the  discrepancy,  and  since  all 
the  other  data  follow  a  predictable  trend,  it  is  considered  an 
outlier . 
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Figure  26  Compactness  Factors  for  All  SUs 


The  entire  set  of  compactness  factor  data  is  plotted  in 
Fig.  28  ,  along  with  a  relationship  for  determining  the 

compactness  factor  as  a  function  of  the  reduction  factor. 
Although  the  spread  in  the  data  is  ±  50.5%,  the  trend  is 
for  smaller  compactness  factors  whenever  rolling  is  done. 

In  an  effort  to  reduce  the  spread,  the  reduced 
compactness  factor  was  plotted  against  the  reduction  factor  in 
Fig.  29.  The  reduced  compactness  factor  is  defined  as  the 
ratio  of  the  compactness  factor  to  the  compactness  factor  for 
unrolled  screens  based  on  empirical  relations  for  Nusselt 
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Compactness  Factor  for  Each  Mesh  Size 


■  100MESH- 

10-^1 - 1 - ^ ^ - 1 - 1 - 1 - ' - 

0.4  0.5  0.6  0.7  0.8  0.9  1  1.1 

Reduction  Factor  -  RF 

Figure  27.  Compactness  Factors  for  Each  Screen  Size 

number  and  friction  factor.  The  best  fit  for  the  data  is  the 
curve,  Reduced  CF  =  Al  Re®^,  where  in  this  case  Al  =  1.22 
±0.061  and  Bl  =  1.65  ±0.211.  The  spread  is  still  moderate  (a 
standard  deviation  of  21.5%  from  the  best  curve  fit),  but  the 
trend  toward  reduced  compactness  factors  as  the  screen 
thickness  gets  smaller  is  still  clear.  This  is  the  most 
important  result  of  the  current  research.  Other  results  are 
also  interesting. 
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Figxire  28.  Global  Compactness  Factor  vs  Reduction 

Factor 


Criteria  Comparisons 

Some  statements  are  made  in  Chapter  III  about  the 
relative  merits  of  the  three  criteria  used  to  select  the  best 
heat  transfer  coefficient.  The  objective  of  this  section  is 
to  show  the  experimental  evidence  which  supports  those 
statements . 

To  begin,  Fig.  30  shows  the  outlet  temperatures  from 
experimental  data  as  well  as  from  the  numerical  model  for 
three  heat  transfer  coefficients,  one  based  on  matching  delay 
time  (DT) ,  one  based  on  minimum  root-mean- squared  difference 
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Figure  29.  Reduced  Compactness  Factor  vs  Reduction 

Factor 

(RMSD) ,  and  one  based  on  matching  the  experimental  maximum 
slope  (MS) . 

In  this  case,  the  delay  time  and  root -mean- squared 
difference  results  are  very  close  to  each  other,  and  match  the 
experimental  trace  for  most  of  its  rise  time.  The  late-time 
effects  of  the  tube  cause  the  traces  to  diverge  there.  The 
maximum  slope  results  do  not  match  well  anywhere.  The  kind  of 
results  shown  in  Fig.  30  did  not  happen  every  time,  but  they 
exemplify  the  types  of  problems  one  can  encounter  unless  the 
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improvements  to  the  classical  technique  made  in  the  current 
research  are  included  in  the  approach. 

To  broaden  the  conclusions  of  the  previous  paragraph, 
other  figures  are  given  below  which  compare  the  delay  time, 
minimum  root-mean-squared  difference,  and  maximum  slope 
approaches . 

The  first  comparison  is  shown  in  Fig.  31  where  the  delay 
time  and  root-mean-squared  difference  results  are  plotted 
against  one  another.  In  general,  they  match  each  other  well 
with  a  mean  difference  of  10.4%.  They  should  match  each  other 
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well  if  the  numerical  model  is  correct  during  the  time  when 
the  root-mean-squared  differencing  is  accomplished.  The 
problem  becomes  choosing  the  appropriate  time  frame.  On  the 
first  run  through  the  data  collected  for  this  research,  the 
root-mean- squared  differencing  was  done  for  the  whole  period 
between  the  beginning  of  the  rise  in  the  outlet  temperature 
trace,  ^2'  to  the  end  of  the  transient,  Tf.  Many  of  the  root- 
mean-squared  difference  heat  transfer  coefficients  were  badly 
in  error  because  the  experimental  outlet  temperature  traces 
diverged  from  the  current  numerical  model  results  at  late 
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times  (as  shown  above)  .  Only  when  the  root-mean-squared 
differencing  time  frame  was  reduced  to  the  period  between  the 
end  of  the  delay  time  to  one  matrix  thermal  time  constant 
after  the  occurrence  of  the  maximum  slope,  did  the  root-mean- 
squared  difference  results  match  the  delay  time  results 
closely.  This  time  frame  is  closer  to  when  the  physics  and 
the  model  match  each  other. 

The  next  comparison  is  made  between  the  delay  time 
results  and  those  for  the  theoretical  maximum  slope  in 
Fig.  32.  As  expected,  the  maximum  slope  approach 

underestimates  the  correct  answer  for  the  heat  transfer 
coefficient  for  the  reasons  given  in  Chapter  III.  The 
literature  shows  a  general  problem  in  accurately  matching 
regenerator  heat  transfer  in  numerical  models  of  regenerative 
systems  (Barnes,  1986:513;  Hamaguchi,  et  al,  1991;  Hutchinson, 
1987;  Tew,  1988) .  The  accepted  conclusion  is  that  oscillating 
flow  causes  an  increase  in  heat  transfer  coefficient,  but  some 
authors  who  have  studied  the  problem  conclude  that  the  heat 
transfer  coefficient  is  not  a  function  of  oscillation 
frequency  (Koester,  et  al  1990) .  Perhaps  the  heat  transfer 
coefficient  data  which  is  being  used  in  current  models  reflect 
the  underestimation  of  the  maximum  slope  method. 

The  third  comparison  shown  in  Fig.  3  3  is  between  the  heat 
transfer  coefficient  based  on  delay  time  and  that  based  on  the 
maximum  slope  of  the  outlet  temperature  trace  calculated  by 
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the  current  numerical  model.  This  approach  is  different  from 
the  analytical  approach  due  to  Pucci,  et  al .  It  takes 
advantage  of  the  contention  that  the  maximum  slope  of  the 
outlet  temperature  trace  is  a  unique  function  of  the  heat 
transfer  coefficient.  The  experimentally  measured  maximum 
slope  is  used  to  choose  the  heat  transfer  coefficient  by 
requiring  the  maximum  slope  of  an  outlet  temperature  trace 
calculated  by  tsie  to  match  the  experimental  one.  Compared  to 
the  delay  time  results  for  Nusselt  number,  eleven  of  the 
twenty  eight  test  runs  match  within  50%.  Two  comparisons  are 
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100%  high,  and  eight  of  the  results  are  66%  lower  than  the 
delay  time  results.  The  balance  of  the  seven  test  runs  did 
not  converge.  A  value  for  the  heat  transfer  coefficient  could 
always  be  calculated  with  the  delay  time  criteria.  Both  the 
large  scatter  and  the  inconsistency  of  using  a  maximum  slope 
make  it  an  inferior  criteria  for  calculating  the  heat  transfer 
coefficient . 

The  final  data  comparison  is  shown  in  Fig.  34.  It 
compares  the  value  of  the  heat  transfer  coefficient  for  the 
analytical  maximum  slope  case  and  the  numerical  maximum  slope 
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case.  In  all  cases,  the  analytical  maximum  slope  predicts  a 
smaller  value  for  the  heat  transfer  coefficient  than  one 
determined  by  an  experimentally  measured  maximum  slope.  This 
was  expected  for  the  reasons  given  in  Chapter  III. 

Effectiveness 

In  Fig.  35,  the  effectiveness  (Eq  (26)),  and  in  Fig.  36, 
the  reduced  effectiveness,  is  shown  as  a  function  of  the 
reduction  factor.  Effectiveness  is  the  dominant  parameter  in 
cryocooler  systems  performance  as  mentioned  in  Chapter  III. 
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Figure  35.  Effectiveness  vs  Reduction  Factor 


The  figures  show  a  reduction  in  effectiveness  whenever 
the  reduction  factor  is  decreased.  In  Fig.  35  the  best  curve 
fit  is  SPEFF  =  A2  Re®^,  where  A2  =  0.89  ±0.018  and  B2  =  0.20 
±  0.059.  The  effect  is  small  at  reduction  factors  of  0.7  and 
0.85,  and  falls  to  about  a  13%  decrement  at  a  reduction  factor 
of  0.5.  Again,  these  results  would  rule  out  any  benefit  of 
reducing  the  screen  thickness  by  50%.  The  interesting  thing 
about  the  figures  is  that  one  would  expect  the  effectiveness 
to  increase  by  making  pore  sizes  smaller  because  of  the  larger 
surface  area  to  volume  ratio.  Rolling  the  screens  does  not 
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change  the  material  properties  of  the  gas  or  matrix,  hence 
another  reason  must  exist  for  the  decrement.  Rolling  the 
screens  does  reduce  the  surface  area  of  the  matrix  since 
abutted  round  wires  have  small  areas  of  contact  while  flat 
surfaces  cover  flat  surfaces  completely.  In  Chapter  II,  large 
surf ace-area-to-volume  ratios  were  credited  with  the  large 
effectiveness  of  wire-screen  matrices,  and  the  area  is 
diminished  when  flat  surfaces  are  introduced  by  rolling  the 
screens . 
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Special  SU 


In  Chapter  III,  SU  #13  was  described.  It  differs  from 
the  rest  of  the  SUs  since  it  is  made  from  a  combination  of 
different  sizes  and  reduction  factors  for  the  screens.  Two 
test  runs  were  accomplished  with  this  regenerator  with  results 
shown  in  Table  4 . 


Table  4.  Comparison  of  S/U  #13  to  the  Unrolled  Case 


Test  Run 

25 

26 

Experimental  / 

1.304 

1.162 

Calculated  / 
(Unrolled  Screens) 

0.927 

0.820 

Per  Cent  Difference 

+  40.7% 

+41.7% 

Experimental  SPEFF 

0.668 

0 . 612 

SPEFF  (Unrolled) 

0.782 

0.725 

Per  Cent  Difference 

-  14.6% 

-  15.6% 

The  purpose  of  these  runs  was  to  gather  data  for  a  truly 
nonhomogeneous  matrix  and  compare  it  with  unrolled  screens. 
Only  data  for  the  friction  factor  is  given  since  the  numerical 
model  would  have  to  rely  on  average  values  for  important 
parameters  to  calculate  the  heat  transfer  coefficient  and  the 
results  may  be  unreliable.  In  both  test  runs,  the  results  for 
the  friction  are  about  40%  higher.  Since  there  is  no 
mechanism  for  enhancing  the  heat  transfer,  the  compactness 
factor  for  this  nonhomogeneous  matrix  must  be  much  smaller 
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than  that  for  a  matrix  with  unrolled  screens  and  the  same 


porosity,  specific  surface  area,  and  mass  flow  rate. 
Additionally,  the  SPEFF  (a  measure  of  effectiveness)  is  shown 
to  be  about  15%  lower  for  both  cases.  While  this  study  is  not 
exhaustive,  it  gives  further  testimony  that  introducing 
heterogeneity  to  the  matrix  causes  a  decrease  in  the 
compactness  factor. 

Reynolds  Analogy 

Due  to  the  lack  of  heat  transfer  data  for  regenerators, 
some  researchers  resort  to  measuring  the  friction  factor  and 
applying  Reynolds  Analogy  (Radebaugh  and  Louie,  1986:180)  to 
determine  the  heat  transfer  coefficient.  For  Reynolds  Analogy 
to  apply  to  flow  in  a  regenerator,  the  compactness  factor 
should  be  a  constant  (Incropera  and  Dewitt,  1985:  272-275) . 
If  the  results  for  flow  in  a  tube  are  used,  Reynolds  analogy 
states  the  compactness  factor  should  be  0.125  (Holman,  1968: 
146)  .  For  the  unrolled  screens,  this  is  not  a  bad  guess 
(16.7%  high  versus  the  average  for  the  unrolled  values)  as 
shown  in  Fig.  37.  But  the  analogy  breaks  down  at  smaller 
reduction  factors,  reaching  a  mean  value  of  0.031  for 
RF  =  0.5.  Reynolds  analogy  works  in  flow  regimes  where  either 
surface  effects  dominate  or  where  inertial  effects  dominate 
(ibid:  274) .  But  in  a  regenerator  where  both  are  important, 
the  analogy  predicts  the  wrong  trend. 
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Figure  37 .  Reynolds  Analogy  Revisited 


Discussion 

In  this  discussion  of  the  results,  three  topics  are 
addressed.  First,  the  uncertainties  for  the  measurements  and 
parameters  is  reported.  Second,  how  changes  in  the  most 
uncertain  parameters  impacts  the  results  is  examined.  The 
section  ends  with  a  summary  of  the  results . 

Uncertainties .  The  uncertainties  for  important 
parameters  and  measurements  made  during  the  current  research 
are  shown  in  Table  5,  along  with  the  source  of  the  estimated 
value.  The  calibrations  for  the  mass  flow  rate,  temperature 
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measurements,  and  pressure  measurements  are  discussed  in 
Appendix  D.  The  three  worst  uncertainties  are  for  the  density 
of  the  steel,  the  mass  flow  rate,  and  the  surface  area  of  the 

,  O 

matrix. 


Table  5.  Uncertainties  in  Measurements 


Quantity 

Uncertainty 

Source 

Specific  Heats,  Cp/Cj^ 

±2.5% 

literature 

Viscosity,  |l 

±3.3 

II 

Porosity,  O 

±0.03 

worst  case 

Reduction  Factor 

±0.0014 

II 

Length,  EL 

±1.3% 

II 

Area  of  Tube,  A 

±1.7% 

Measured 

Pore  Size,  PS 

±1.3  )im 

II 

Wire  Diameter, 

It 

n 

Matrix  Mass,  W^u 

±0.001g 

II 

Delay  Time,  SEDT 

±0 . OOlsec 

II 

Matrix  Density, 

±  10% 

lit. /self  test 

Gas  Density,  p 

±  1% 

calculated 

Mass  Flow  Rate,  m 

+7.5%,  -0% 

calibration 

Matrix  Surface 
Area , 

+0,  -20% 

calculated 

Temperature,  T 

1+ 

o 

to 

0 

n 

calibration 

Pressure,  P 

±1%  of  full  scale 

calibration 

Effects  of  Uncertainties.  The  effect  on  the  data  for 
each  of  the  three  worst  uncertainties  is  discussed  in  this 
section.  The  first  is  the  uncertainty  of  the  density  of  the 
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steel.  Any  error  in  this  material  property  would  affect  all 
the  data  runs  similarly,  hence,  the  any  conclusions  would  not 
change  since  they  are  based  on  a  comparison  between  the 
characteristics  of  rolled  and  unrolled  SUs . 

The  measured  mass  flow  rate  uncertainty  was  calibrated  to 
be  ±  7.5%  for  the  lowest  Reynolds  number  and  ±  4.0%  at  the 
highest.  If  the  flow  velocity  is  really  larger  than  measured, 
the  effect  on  the  compactness  factor  can  be  seen  by  studying 
the  Armour  and  Cannon  (1968)  relation  for  friction  and  the 
curve-fit  for  the  Nusselt  number  given  above.  Since  the 
Reynolds  numbers  would  be  larger,  the  friction  factor  would  be 
smaller.  The  Nusselt  number  would  be  larger  proportional  to 
ReO.62,  lienee  the  Stanton  number  would  be  smaller  proportional 
to  Re'°'^®.  The  ratio  between  the  Colburn  factor  and  the 
friction  factor  would  make  the  uncertainty  in  the  compactness 
factor  ±  7.8%  for  the  lowest  flow  rates  and  ±  2.0%  for  the 
largest  flow  rates .  Since  these  spreads  are  well  below  the 
spread  in  the  experimental  data  (±  37.4%),  the  results  are 
unchanged . 

The  largest  uncertainty  is  in  the  surface  area  of  the 
matrix.  This  problem  was  addressed  by  first  assuming  there 
was  no  surface  area  lost  due  to  rolling.  The  effect  on  the 
Nusselt  number  was  to  make  it  smaller  by  a  mean  difference  of 
6.1%,  compared  to  the  subsequent  analysis  which  included  the 
area  correction.  The  friction  factor  was  unaffected.  The 
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treatment  of  the  area  was  somewhat  conservative  since  it  was 


assumed  all  flattened  areas  were  no  longer  available  for  heat 
transfer  or  surface  friction  exchanges.  Hence,  the  final 
answer  for  the  compactness  factor  lies  somewhere  between  the 
6.1%  difference  between  the  calculated  Nusselt  numbers  with 
and  without  area  change,  which  is  within  the  ±  37.4%  spread  in 
the  plot  of  compactness  factor  vs  reduction  factor. 

Summary .  The  results  presented  in  this  chapter  give  a 
convincing  explanation  for  the  conclusion  that  flattened 
screens  causes  the  ratio  of  heat  transfer  coefficient  to 
friction  factor  to  decrease.  Confidence  in  the  results  is 
warranted  because  of  the  favorable  comparisons  to  other  data, 
the  repeatability  of  the  test  runs,  the  improvements  made  in 
the  methodology,  the  small  uncertainties  in  the  measurements, 
and  the  immutability  of  the  results  given  a  worse  case  change 
in  measurements . 

The  data  show  that  the  heat  transfer  coefficient  is 
unaffected  by  the  rolling  of  the  screens  when  the  reduction  in 
the  wetted  surface  area  caused  by  rolling  the  screens  is  taken 
into  account.  No  enhancement  to  the  heat  transfer 
characteristics  of  the  regenerator  is  created  by  introducing 
flattened  surface  to  the  geometry. 

The  friction  factor,  however,  appears  to  increase 
marginally  for  the  larger  reduction  factors,  but  substantially 
for  the  50%  reduction  case.  As  mentioned  in  Chapter  II,  there 
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are  two  components  to  the  pressure  drop,  a  surface  shear  and 
an  inertial  component.  Although,  measurements  of  these  two 
components  were  not  made,  an  argument  for  the  observed 
pressure  drop  behavior  can  be  made.  Reducing  the  surface  area 
would  reduce  the  surface  shear,  but  decreasing  the  size  of  the 
pores  and  introducing  sharp  edges  into  the  geometry  would 
increase  the  inertial  component  of  the  pressure  drop.  At  the 
higher  reduction  factors,  RF=0.7  or  especially  RF  =  0.85,  the 
inertial  component  of  the  pressure  drop  increases  to 
compensate  for  the  reduction  in  surface  drag  forces  due  to  a 
decrease  in  surface  area  which  leaves  the  total  pressure  drop 
nearly  constant.  At  a  reduction  factor  of  0.5,  there  is  a 
large  increase  in  inertial  pressure  drop.  The  pore  size 
becomes  so  small,  and  the  velocity  increases  so  much,  the 
onset  of  compressibility  effects  occurs  which  causes  the 
friction  factor  to  grow  much  larger. 

The  combination  of  the  results  above  leads  to  the 
conclusion  the  compactness  factor  decreases  whenever  rolling 
the  screens  is  done.  Hence,  the  hypothesis  presented  in 
Chapter  I  is  proven:  Using  pressed  screens  to  construct  the 
regenerator  matrix  decreases  the  volume  of  the  regenerator, 
but  by  doing  so,  the  compactness  factor  is  also  reduced. 
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VI .  Conclusions  and  Recommendations 


The  purpose  of  this  last  chapter  is  to  suininarize  the 
current  research,  and  make  suggestions  for  further 
investigations . 

The  objective  of  the  research  was  to  measure  the  heat 
transfer  and  friction  factor  characteristics  for  a  typical 
stacked-screen  regenerator  which  has  had  its  geometry  altered 
by  reducing  the  thickness  of  its  screens.  Regenerative  heat 
exchanger  technology  is  important  for  space-based  energy 
conversion  systems,  and  any  increase  in  performance  which  can 
be  achieved  by  incorporating  easily  achieved  design  features 
should  be  put  to  full  advantage.  Some  researchers  have 
reported  an  improvement  in  system-level  performance  of  an 
engine  by  rolling  the  screens  in  a  stacked-screen  matrix. 
Although  this  certainly  reduced  the  dead  volume,  they  did  not 
report  whether  the  improvement  in  system- level  performance  was 
due  to  the  decrease  in  dead  volume,  or  to  an  improvement  in 
the  operating  characteristics  of  the  regenerator,  represented 
by  the  compactness  factor,  or  to  a  combination  of  both.  Since 
the  benefits  of  reduced  dead  volume  have  already  been 
documented,  a  study  of  the  effect  on  the  compactness  factor 
was  undertaken  for  a  range  of  geometries  and  flow  conditions 
found  in  current  cryocooler  systems. 

Rolling  the  screens  did  reduce  the  volume  of  the 
regenerator,  but  the  results  of  the  current  research  show 
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rolling  the  screens  did  nothing  to  enhance  the  heat  transfer. 
However,  the  pressure  drop  increases,  particularly  for 
reduction  factors  of  50%.  At  larger  reduction  factors,  around 
0.85,  the  decrease  in  heat  transfer  and  increase  in  pressure 
drop  are  small,  which  left  the  compactness  factor  nearly 
unchanged.  But  the  effectiveness  data  indicate  any  reduction 
in  screen  thickness  would  degrade  the  effectiveness,  the 
dominant  factor  for  determining  cryocooler  coefficient  of 
performance.  With  the  friction  and  heat  transfer  data 
presented  in  Chapter  IV,  a  cryocooler  designer  could  do  a 
trade-off  study  which  includes  the  influence  of  effectiveness, 
pressure  drop,  and  volume  reduction,  and  determine  which 
combination  of  flow  conditions,  geometry,  and  reduction  factor 
gives  the  optimum  performance. 

The  crux  of  the  matter  with  regard  to  the  compactness 
factor  is  how  the  geometry  of  the  matrix  changed  after  the 
screens  had  been  rolled.  Rolling  the  screens  caused  the 
smooth,  round  surfaces  of  wires  to  be  replaced  with  flattened 
areas  with  abrupt  edges.  The  flattened  areas  abut  one  another 
causing  the  available  wetted  surface  area  to  decrease. 
Reduced  wetted  surface  area  reduces  the  total  heat  transfer. 
But  the  smaller  pore  sizes  and  sharp  edges  causes  the  total 
pressure  drop  to  increase.  The  overall  effect  of  the  geometry 
change  is  to  decrease  the  compactness  factor. 
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The  conclusions  concerning  heat  transfer  are  due  in  part 
to  the  improved  method  for  determining  the  heat  transfer 
coefficient  in  a  porous  medium  used  in  the  current  research. 
The  measured  inlet  temperature  profile  is  used.  The  sponge 
effect  delay  time  during  which  the  numerical  model  and  the 
physics  of  the  flow  more  closely  match  one  another,  is  a 
characteristic  of  the  matrix  and  is  the  criteria  for  choosing 
the  heat  transfer  coefficient  .  Also,  the  significant  effect 
of  the  tube  surrounding  the  matrix  is  included  in  the  model. 
These  modifications  in  the  method  resulted  in  generally  higher 
values  for  the  Nusselt  number  than  occurs  in  other  transient 
step-change  techniques . 

Further  study  in  four  areas  warrants  attention.  First, 
the  same  apparatus,  data  reduction,  and  numerical  model  could 
be  used  with  an  oscillating  flow  of  gas  to  determine  if 
oscillation  frequency  is  a  factor  in  the  heat  transfer  or 
pressure  drop  characteristics.  This  would  be  particularly 
interesting  because  the  results  could  be  obtained  with  the 
same  apparatus  and  the  same  technique,  a  condition  not  found 
in  current  experimentation.  Second,  since  there  are  ways  to 
increase  the  surface  area  of  the  wire  before  it  is  rolled, 
e.g.,  by  etching  or  drawing  small  grooves  in  the  wire,  the 
effectiveness  of  a  regenerator  might  be  increased,  even  for  a 
rolled  screen  configuration.  Thirdly,  since  predicting  the 
characteristics  of  regenerators  which  are  not  constructed  in 
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a  homogeneous  fashion  is  impossible,  a  further  study  of  sponge 
effect  delay  times  could  determine  the  usefulness  of  comparing 
regenerators  based  solely  on  delay  time  characteristics  as  a 
figure  of  merit.  Finally,  analytical  studies  which  include 
the  compactness  factor  and  Reynolds  analogy  in  lieu  of 
experimental  heat  transfer  data  to  estimate  the  performance  of 
porous  medium  regenerators  could  be  repeated  to  determine  how 
the  decreasing  trend  in  compactness  factor  changes  their 
results.  In  particular,  since  heat  transfer  and  friction 
factor  data  are  now  available,  a  study  which  incorporates  all 
the  influences  of  effectiveness,  pressure  drop,  and  reduced 
dead  volume  on  the  system-level  performance  should  be 
undertaken  to  determine  the  best  configuration  for  a  stacked, 
wire-screen  regenerator  with  reduced  thickness  screens. 
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TEST  PROCEDURE 


(CAPT  TIM  MURPHY, AFIT  DSY9 6M, REGENERATOR  HEAT  TRANSFER  STUDY) 
RIG  SET  UP 


-  TURN  ON  COOLER/HOT  PLATE: 

-  MAKE  SURE  ENOUGH  HE  FOR  TEST: 

-  TURN  ON  DC  POWER  SOURCE _  VI _  IFA  100 _  VALIDYNE. 

-  CHECK  1  VDC  ±  0.001  ON  THERMISTOR  CIRCUIT 

-  SAFETY  EQUIPMENT  IN  PLACE _ 

-  LEAK  CHECK  WITH  SOAPY  WATER  FOR  S/U  SWITCHOUT _ 

-  WAIT  ONE  HOUR  FOR  STEADY  STATE _ 


TEST  RUN 


-  FILL  IN  DATA  SHEET  FOR  THIS  RUN _ 

-  RECORD  ROOM  TEMPERATURE/ PRESSURE  ON  .VI _ 

-  MAKE  SURE  SAMMPLE  FREQUENCY  AND  BUFFER  ARE  SET  TO  1000 /SEC _ 

-  ENSURE  REGEN. RAW  IS  EMPTY  OF  OLD  TEST  RUN  DATA _ 

-  CLOSE  VENT  OUTLET  VALVE;  4 -WAY  VALVE  TO  HOT  FLOW 

-  RUN  HOT  HE  UNTIL  T1/T2  EQUILIBRATE  THEN  SHUT  OFF  THE  VALVE _ 

-  SET  TANKS  BACK  PRESSURE  TO  APPROPRIATE  _  psia 

-  ZERO  OUT  DPVEN _  DPSU _  EPl _ PRESSURE  TRANSDUCERS 

-  PUT  4 -WAY  VALVE  TO  COLD  FLOW  AND  OPEN  UPSTREAM  VALVE _ 

-  ADJUST  BOTH  VALVE  UNTIL  PI  AND  MDOT  ARE  ABOUT  RIGHT _ 

-  PUT  4 -WAY  TO  HOT  FLOW  AND  MATCH  COLD  FLOW  PI  AND  MASS  FLOW _ 

-  TURN  4 -WAY  BACK  TO  COLD  FLOW  AND  WAIT  FOR  TEMPERATURE  TO 

EQUILIBRATE  AT  ABOUT  10°C _ 

-  LET  HOT  HE  TRICKLE  UNTIL  TEMPERATURE  IN  VENT  IS  APPROXIMATELY 

3  5°C _ 

-  HIT  "SEND  DATA  TO  FILE"  ON  .VI _ 

-  WAIT  THREE  SECONDS;  TURN  4 -WAY  VALVE  TO  HOT  FLOW _ 

-  WAIT  20  SECS;  NOTE  PRETRANSIENT  MEASUREMENTS _ 

-  TURN  OF  GAS  FLOW  HOT _  COLD _ 

-  NOTE  POST  TRANSIENT  MEASUREMENTS _ 

-  OBSERVE  DATA  TRACE  FOR  QUALITY;  REPEAT  IF  NECESSARY _ 

-  RENAME  REGEN. RAW  FOR  THIS  TEST  (t**r**.raw) _ 

-  MAKE  NOTES  ABOUT  THE  TEST  RUN  ON  THE  DATA  SHEET _ 

-  ENTER  TEST  RUN  IN  THE  LOG  BOOK _ 


DATA  SHEETS 


DATE: _  TIME: _ 

TEST  #: _  S/U: _  PORE 

SIZE: _ 

REps  : _  MASS  FLOW  RATE: _ 

ROOM  BAROMETRIC  PRES  (psia)  : _  ROOM  TEMP,  (deg 

C)  : _ 

HE  TANK  PRESSURE  (psia)  : _  BACK 

PRESSURE  (psia)  :_ _ 

DATA  FILE  NAME: _ 

ANEMOMETER  OSCILLATIONS  ?: _ 

PRETRANSIENT:  DPSU _  DOTM _  DELTA 

T1/T2 _ 

PI _ 

POST-TRANSIENT:  DPVEN _  DPSU _ 

PI _ 


COMMENTS : 
OK? 


DELTA  T1/T2 


TEST  QUALITY 
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function  range=trange (tlt2, t) 

%The  purpose  of  this  function  is  to  take  in  a  range  of  independent  variables,  tr 
%and  a  function,  temp  ,of  the  independent  variable  and  determine  when  a  'change 
%is  occuring,  e.g.  when  the  slope  of  the  function  is  more  than  1%  different 
%  from  the  mean  slope  for  the  range.  This  will  define  a  range  of  the 
%  independent  variable  for  which  a  transient  is  defined. 

g, 

0 

delt=diff (t (1 :2) ) ; 
tl=tlt2  ( :, 1) ; 
t2=tlt2  (:,2) ; 
b=firl (100, .006) ; 
tlf=filtfilt (b, 1, tl) ; 
t2f=filtfilt (b, l,t2) ; 
l=length  (tlf ) ; 
lder=l-l; 
lz=lder-220; 
lzl=lder-121; 
derl=diff (tlf) ./diff (t)  ; 
der2=diff (t2f) ./diff (t) ; 
tav=mean (tlf (21 : 121) ) ; 
derav=mean (der2 (lz;l2l) ) 

%plot (der2) 

% 

ti=zeros (1, 1) ; 
for  i=l:l 

if  (tlf(i)-tav)  >  1.0; 
ti(i)=-1.0/t(i) ; 
end 

ti  (i) =ti  (i) +50 . ; 
end 

[time  k]=sort(ti); 
indi=k (1) ; 
tri=t (indi) ; 

% 

tf=zeros  (1, 1) ; 
for  j=lz:-l;indi 

if  der2(j)  >  5.0; 
tf (j)=-1.0*t(j); 
end 

end 

[time  k]=sort(tf); 
indf=k(l) ; 
trf=t(indf); 
range=[indi  indf ] ; 

plot  ( (tri  :delt  :trf ) ,  tl  (indi : indf ) ,  (trirdelt  :trf ) ,  t2  (indi : indf) ) 


function  inslope=inaxslope  (p,  tr) 

%This  function  takes  in  the  coefficients  of  a  polynomial, p,  and  a  range  of 

%independent  variables, tr,  and  returns  the  value  of  the  maximum  slope  in  that  range 

%as  well  as  the  location. 

ord=length (p) -1; 

dor=ord-l 

pl=[p(l:ord) ] ; 

n=  (ord:-l ; 1) ; 

pd=pl . *n; 

der=polyval (pd,  tr) ; 
lder=length (der) ; 
plot (der) 

[md  loc]=sort (der' ) ; 
lms=tr (loc (Ider) ) ; 
ms=md(lder) ; 
mslope=[ms  1ms] 


function  bestorder=bestpoly (temp,  tr) 

%This  function  takes  vectors  of  dependent  variables,  teriro,  and 
%independent  variables, tr,  and  returns  the  best  order  of  the 
%polynomial  which  fits  the  temp  data  in  the  range  tr  based  on 
%the  minimum  value  of  the  root  mean  squared  difference, 
for  i=l:27 

p=polyfit  (tr,  ten^D,  i) ; 
f=polyval (p,  tr) ; 
dif=teii5i-f  ; 
dif=dif ' ; 
sdif=dif .  ''2; 
cdif=sum(sdif) ; 
msdif=cdif/length(dif) ; 
rms=sqrt (msdif ) ; 

M=mean (dif ) ; 

S=std(dif) ; 

MS  (i, l)=rms; 

MS(i,2)=S; 

MS (i, 3)=M; 
end 

ams=abs (MS) 

[ms  j]=sort (ams) ; 
ms; 

check= j (1, 1) - j  (1,2); 

3 

minrmsdif=ms (1,1) 
order= j (1,1) 

bestorder=polyfit  (tr, temp, order) 


function  dt=t2cit  (t2,  T,  II,  IF) 

%This  m.file  takes  in  a  data  file  containing  the  temperatures  and  times 
%of  a  thermal  transient  and  returns  dt,  the  time  it  takes  for  the 
%temperature  to  change  "significantly"  from  its  initial  value. 

*0 

b=firl(20, .006) ; 
t2f=filtfilt (b, 1, t2) ; 
l=length (t2) ; 

11=1-100; 
ti=ones  (1, 1) ; 

9- 

o 

for  i=II:IF 

tav=mean (t2f (II : i) )  ; 
if  t2f (i+50) -tav  >  0.4; 
ti (i) =-l . 0/i; 
end 

end 

Q, 

0 

[time  k]=sort(ti); 
i2=k{l); 
tau2=T (i2) ; 
dt=[i2  tau2] 

plot  (T  (i2  :  IF) ,  t2f  (i2  :  IF) ) 


function  teittp=filtlt2  (tl,  t2) 

% 

%This  m.file  takes  in  the  values  of  the  inlet  temperature, tl, 

%and  the  outlet  temperature,  t2,  and  filters  them  with  a  100 

%order  fir  filter,  tl  and  t2  are  two  vectors  containing  experimenetal 

%data . 

*0 

b=firl(20, .006) ; 
tlfir20=filtfilt  (b, 1, tl) ; 
t2fir20=filtfilt  (b, 1, t2) ; 
temp= [tlfir20  t2fir20]; 
plot (tlfir20) 
hold 

plot (t2fir20) 
hold  off 
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PROGRAM  TSIE 

C  THIS  IS  A  COMPUTER  PROGRAM  TO  DETERMINE  THE  TRANSIENT  GAS 

C  TEMPERATURE  FOR  A  WIRE-SCREEN  REGENERATOR  SUBJECTED 

C  TO  A  QUASI  STEP  CHANGE  IN  TEMPERATURE  AT  THE  INLET. 

C  IT  USES  A  STRICTLY  EXPLICIT  FINITE  DIFFERENCE 

C  SCHEME  AT  EACH  TIME  STEP  AND  INCLUDES  THE  INTERACTION  BETWEEN 

C  THE  GAS  AND  THE  TUBE  ALTHOUGH  NOT  CONDUCTION  BETWEEN  THE  TUBE 

C  AND  THE  MATRIX. 

C  THE  WORKING  FLUID  IS  HELIUM  AT  ROOM  TEMPERATURE. 

C  THE  MASS  FLOW  RATE  IS  GIVEN  AND  ASSUMED  STEADY. 

C  THE  ENERGY  EQUATION  FOR  INCOMPRESSIBLE  FLOW  IS  USED  TO  CALCULATE 

C  THE  TEMPERATURE  OF  THE  EXIT  GAS. 

C  DOTM  IS  THE  MASS  FLOW  RATE  (FUNCTION  OF  TIME  ONLY)  . 

C  TGO(I)  AND  TGL(I)  ARE  THE  EXPERIMENTALLY  MEASURED  GAS  TEMPS 

C  AT  THE  INLET  AND  EXIT  OF  THE  REGENERATOR,  RESPECTIVELY . 

C  D  IS  THE  PIPE  DIAMETER;  EL  IS  THE  REGENERATOR  LENGTH. 

C  POR  IS  THE  PORISITY;  I,N  ARE  THE  TIME  AND  SPACE  COUNTERS. 

C  DH  IS  THE  HYDRAULIC  DIAMETER  OF  THE  REGENERATOR  MESH. 

C  DELT,  DELX  ARE  THE  TIME  AND  SPACE  INCREMENTS  OF  THE  MODEL. 

C  SAMPT  IS  THE  EXPERIMENTAL  DATA  SAMPLING  TIME  INCREMENT. 

C  I,J,  AND  N  ARE  THE  TIME  AND  SPACE  COUNTERS. 

C  II  AND  IF  ARE  THE  BEGINNING  AND  ENDING  INDEXES  ASSOCIATED  WITH 

C  TAUO  AND  TAUF  OF  THE  ACTUAL  TRANSIENT. 

C  UMAX  AND  NMAX  ARE  THE  MAXIMUM  TIME  AND  SPACE  INCREMENTS. 

C  ITRAN  IS  THE  NUMBER  OF  TIME  INCREMENTS  DURING  THE  ACTUAL 

C  TRANSIENT,  TAUF-TAUO/DELT . 

C  10  IS  THE  COUNTER  FOR  THE  EXPERIMENTAL  DATA  AT  THE  BEGINNING 

C  OF  THE  CURRENT  I COUNT. 

C  IINT  IS  THE  NUMBER  IF  CALCULATION  TIME  INCREMENTS  BETWEEN 

C  EXPERIMENTAL  DATA  POINTS. 

C  ICOUNT  IS  THE  COUNTER  FOR  THE  CALCULATION  SUBINTERVAL  BErWEEN 

C  DATA  POINTS. 

C  EMU  IS  THE  DYNAMIC  VISCOUSITY;  ENU  THE  KINEMATIC  VISC. 

C  REPS  IS  THE  REYNOLD'S  NUMBER  BASED  ON  PORE  SIZE. 

C  PS  IS  THE  DISTANCE  BETWEEN  WIRES. 

C  DP  IS  THE  WIRE  DIAMETER. 

C  PIT  IS  THE  PITCH  =  DP+PS. 

C  HTC  IS  DEFINED  AS  Q  /  A/  TG-TR 

C  TGC  (N)  ,  TRC  (N)  ,  AND  TTC  (N)  ARE  THE  TEMPERATURES  OF  THE  GAS,  MATRIX 

C  AND  TUBE  WALL  DURING  THE  CURRENT  TIME  STEP  OF  THE  TRANSIENT. 

C  TGP(N),  TRP(N),  AND  TTP (N)  ARE  THE  GAS, MATRIX, AND  TUBE 

C  TEMPERATURES  DURING  THE  PREVIOUS  TIME  STEP  OF  THE  TRANSIENT. 

C  TG2  IS  THE  OUTLET  GAS  TEMPERATURE  DETERMINED  BT  THE  MODEL. 

C  TEMP  IS  THE  INTERPOLATED  INLET  GAS  TEMPERATURE. 

C  HNU  IS  THE  NUSSELT  NUMBER;  NU  =  HTC*DH/KG 

C  R  IS  THE  GAS  CONSTANT  FOR  THE  WORKING  FLUID. 

C  EKG  IS  THE  THERMAL  CONDUCTIVITY  OF  THE  GAS. 

C  EKS  IS  THE  THERMAL  CONDUCTIVITY  OF  THE  TUBE  MATERIAL. 

C  A  IS  THE  TOTAL  PIPE  AREA;  AF  IS  THE  AREA  FOR  FLOW  =  POR*A. 

C  AS  IS  THE  SURFACE  AREA  FOR  HEAT  TRANSFER  PER  UNIT  VOLUME  OF 

C  THE  REGENERATOR  MATRIX. 

C  AST  IS  THE  TUBE  AREA  FOR  HEAT  TRANSFER  IN  ONE  DELX. 

C  ACXT  IS  THE  CROSS  SECTIONAL  AREA  OF  THE  TUBE  WALL. 

C  CP  IS  THE  SPECIFIC  HEAT  OF  THE  GAS  AT  CONSTANT  PRESS. 

C  CV  IS  THE  SPECIFIC  HEAT  OF  THE  GAS  AT  CONSTANT  VOLUME. 

C  GAMMA  IS  THE  RATIO  OF  CP/CV. 

C  SLOPE  IS  THE  MAXIMUM  SLOPE  OF  THE  TEMPERATURE  TIME 

C  CURVE  AT  t  =  ELOC;  LOC  IS  THE  INDEX  FOR  ELOC. 

C  RHO  IS  THE  GAS  DENSITY  =  PM/R/TM 

C  CM  IS  THE  MATRIX  HEAT  CAPACITY;  RHOM  THE  MATRIX  DENSITY. 

C  WSU  IS  THE  MASS  OF  THE  SCREEN  UNIT. 

C  THK  IS  THE  TUBE  THICKNESS. 

C  MESH  IS  THE  MESH  NUMBER  (PER  INCH,  SEE  MIYABE.) 

C  SAMPT  IS  THE  EXPERIMENTAL  DATA  SAMPLING  TIME. 

C  TTOT  IS  THE  TOTAL  DATA  TIME,  i.e.  IMAX  =  TTOT/SAMPT 

C  TAUO  AND  TAUF  ARE  BEGINNING  AND  ENDING  TIMES  FOR  THE  ACTUAL 
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C  TRANSIENT. 

C  TAU2  IS  THE  TIME  FOR  THE  SEDT  TO  END. 

C  T0I,TLI  ARE  THE  INITIAL  TEMPERATURES  AT  THE  INLET  AND  OUTLET. 

C  THETA  IS  THE  ANGLE  OF  CONTACT  BETWEEN  STRANDS  OF  THE  SCREEN. 

C  RF  IS  THE  REDUCTION  FACTOR  FOR  THE  SCREEN,  RF=EL  (act) /ELO  . 

C 

DOUBLE  PRECISION  TGC (155) , TRC (155) , TTC (155)  ,  TGP ( 155) , TRP (155) , 
$TTP (155)  ,DOTM, REPS, HTC, HNU, TGL (20001) , TGO (20001) , TG2 (20001) 
INTEGER  I, N, NMAX, JMAX, LOC, II, IF, 10, IINT, ICOUNT, ITRAN, J 
DOUBLE  PRECISION  R,  RHO,  TA,  TB,  TC,  TD,  CM,  RHOM,  A,  CP,  CV,  GAM,  TIME,  THK, 
$DELT,  DELX,  DH,  AF,  ENU,  EMU,  POR,  EL,  D,  PI,  DP,  PS,  EKG,  EKS,  PIT,  THETA, 

$TE,  TF,  TH,  T  J,  PAV,  TOI,  TLI,  TAV,  SLOPE,  AS,  C,  TTOT,  ELOC,  TAUO,  TAUF, 

$ACXT,  AST,  SAMPT,  TEMP,  WSU,  SLOPEl,MESH,  RF,  TG2AV,  SEDT,  TAU2 
C 

INTRINSIC  COS, MAX, SQRT, ABS,  NINT,  ATAN 

OPEN  THE  DATA  OUTPUT  FILE 
OPEN (50,  FILE='  tss . in' , STATUS=' UNKNOWN' ) 

OPEN(52,FILE=' tss .out' , STATUS=' UNKNOWN'  ) 

OPEN (53, FILE=' tsie.dat', STATUS=' UNKNOWN'  ) 

INPUT  FLOW  PARAMETERS  AND  NUMERICAL  VARIABLES. 

INPUT  THE  WIRE  DIAMETER  AND  PORE  SIZE 
READ(50,102)  DP, PS 
WRITE (53,*)  'DP  =' ,  DP,  'PS  =  '  ,  PS 

INPUT  REGENERATOR  LENGTH  AND  MASS 
READ (50, 102)  EL, WSU 

WRITE  (53,*)  'LENGTH  =  ',  EL,  'WEIGHT  =  '  ,  WSU 

WRITE  (*,*)  '  INPUT  THE  HEAT  TRANSFER  COEFFICIENT, HTC, AND  DELT  ' 
READ(*,*)  HTC,  DELT 

WRITE (53,*)  '  HTC  EQUALS  ',HTC, '  DELT  =  ',DELT 

INPUT  THE  MESH  NUMBER 
READ (50, 103)  MESH 
WRITE (53,*)  MESH,  'MESH' 

THE  TUBE  DIAMETER  AND  THICKNESS  ARE  CONSTANT 
0=0.014859 
THK=0. 000508 

ENTER  THE  NUMBER  OF  SPACE  INCREMENTS 
WRITE (*,*)  'INPUT  THE  NUMBER  OF  SPACE  INCREMENTS' 

READ(*,*)  NMAX 
WRITE (53,*)  '  NMAX  =  ',  NMAX 

R  IS  THE  GAS  CONSTANT  FOR  HELIUM  IN  J/Kg*K 
R=2077. 


DETERMINE  THE  MAXIMUM  THE  SPACE  INCREMENTS. 

DELX=EL/ (NMAX-1) 

ADD  TWO  SPACE  INCREMENTS  TO  THE  TOTAL. 

NMAX=NMAX+2 

EKG  IS  THE  THERMAL  CONDUCTIVITY  OF  HELIUM, W/m-K.  (SEE  KAYS) 
EKG=0.155 

EKS  IS  THE  THERMAL  CONDUCTIVITY  OF  304  SS, W/m-K. (SEE  INCROPERA) 
EKS=14 . 9 

NEED  A  VALUE  FOR  EMU  FOR  THE  HELIUM  (Pa*SEC)  . 

EMU=1.95E-05 


RHO  IS  THE  DENSITY  OF  THE  GAS,  ASSUMED  CONSTANT 
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READ (50, 102)PAV, TAV 

WRITE (53,*)  'PAV  =  ',  PAV,  '  TAV  =  ',  TAV 
RH0= (PAV* 68 95 . ) /R/ (TAV+273 . ) 

DETERMINE  THE  KINEMATIC  VISCOSITY. 

ENU=EMU/RHO 

DEFINE  PI 
PI=4.0*ATAN(1.0) 

NEED  THE  TOTAL  AREA  OF  THE  PIPE. 

A=PI*D**2/4.0 

NEED  THE  CROSS  SECTIONAL  AREA  OF  THE  TUBE  WALL,ACXT. 
ACXT=PI* (D+THK/2 .0) *THK 

NEED  A  VALUE  FOR  THE  THERMAL  CAPACITY  OF  THE  MATRIX 
MATERIAL  IN  J/KG/K 
CM=477 . 0 

NEED  A  VALUE  FOR  THE  DENSITY  OF  THE  MATRIX  IN  KG/M3 . 
RHOM=7100. 

NEED  A  VALUE  FOR  THE  HEAT  CAPACITY  OF  THE  GAS  IN  J/KG/K. 
CP=5193. 

NEED  A  VALUE  FOR  GAMMA  AND  CV. 

CV=CP-R 

GAM=CP/CV 

NEED  A  VALUE  FOR  HEAT  TRANSFER  AREA,  AS.  (SEE  MIYABE /ARMOUR) 
ENTER  THE  REDUCTION  FACTOR. 

READ (50, 103)  RF 
PIT=DP+PS 

DETERMINE  THE  CONTACT  ANGLE. 

THETA=ATAN (DP/PIT) 

WRITE (53,*)  '  RF=  ',RF 

ASP=PI* (SORT (PIT**2+DP**2)-  (DP*THETA**2/PI) )  /PIT**2 

CORRECT  FOR  THE  REDUCTION  FACTOR  (SEE  NOTES.) 

IF  (RF.LT.0.65)  THEN 
ASP=ASP*0.8 
GO  TO  15 

ELSEIF  (RF.LT.0.8)  THEN 
ASP=ASP*0.97 
GP  TO  15 

ELSEIF  (RF.LT.0.95)  THEN 
ASP=ASP *0.9992 
GO  TO  15 

ELSE 

GO  TO  15 

END  IF 
CONTINUE 

AS=ASP*A* (DELX/RF) 

ATOT=AS*EL/DELX 

WRITE (53,*)  'HEAT  TRANSFER  AREA  EQUALS  ',AS,'  m2' 

DETERMINE  THE  POROSITY  (SEE  MIYABE  1982) . 

POR=l . 0- (WSU/RHOM/A/EL) 

WRITE (53,*)  POR, 'POR' 

DETERMINE  THE  FLOW  AREA. 

AF=POR*PI*D**2/4 . 0 

NEED  THE  TUBE  HEAT  TRANSFER  AREA,  AST. 

AST=P I *D  *DELX*POR 


c 

C  DEFINE  THE  HYDRAULIC  DIAMETER. 

DH=POR/ (AS/A/DELX) 

WRITE (53,*)  'THE  HYDRAULIC  RADIUS  IS  ',DH, '  m2' 

C 

C  INPUT  THE  SAMPLING  FREQUENCY  (  =  1/SAMPT) . 

READ (50, 102)  SAMPT, TTOT 

C  WRITE  (53,*)  SAMPT, TTOT, 'SAMPT  and  ttot' 

C  INPUT  THE  TOTAL  TIME  FOR  THE  TRANSIENT  -  IN  SECONDS 

JMAX=NINT (TTOT/SAMPT) 

C  WRITE (53,*)  UMAX,'  JMAX=  ' 

105  FORMAT (19) 

C 

C  DEFINE  THE  COURANT  NUMBER  C. 

C=DELX/DELT 

WRITE (53,*)  '  THE  COURANT  NUMBER  IS  ',C 
C 

N=1 

J=1 

1=1 

C 

C  READ  IN  THE  MASS  FLOW  RATE  AND  THE  INLET  TEMPERATURES  WHICH 

C  APPROXIMATE  A  STEP  CHANGE. 

READ  (50,  103)  DOTM 

WRITE (53,*)  DOTM,  '  =  dotm  ' 

DO  43  J=1,JMAX 

READ (50, 104)  TGO (J) ,TGL(J) 

43  CONTINUE 

C  WRITE(53,*)  TG0(56),TGL(56),'  tgO  (56) ,  tgl (56) ' 

C 

C  DETERMINE  THE  INDEX  OF  THE  BEGINNING  OF  THE  ACTUAL  TRANSIENT,  II, 

C  AND  THE  END  OF  THE  ACTUAL  TRANSIENT,  IF . 

C 

READ (50, 102)  TAU0,TAUF 

WRITE (53,*)  TAU0,TAUF,'  TAU0,TAUF  ' 

II=NINT (TAUO/SAMPT) +1 
IF=NINT (TAUF/ SAMPT) +1 
WRITE (53,*)  II, IF,  '  II, IF  ' 

C 

C  DETERMINE  THE  NUMBER  OF  INCREMENTS  BETWEEN  DATA  POINTS,  IINT,  AND 

C  THE  NUMBER  OF  TIME  INCREMENTS  DURING  THE  ACTUAL  TRANSIENT,  ITRAN. 

IINT=NINT (SAMPT/DELT) 

ITRAN=NINT ( (TAUF-TAUO) /DELT) +1 
C  WRITE (53,*)  IINT, ITRAN, '  IINT, ITRAN  ' 

C 

C  SET  THE  INITIAL  CONDITION  FOR  THE  GAS  AND  MATRIX  AT  A  LINEAR 

C  DISTRIBUTION  BETWEEN  TOI  AND  TLI.  TOI  IS  A  MEAN 

C  INITIAL  INLET  TEMP.  TLI  IS  THE  MEAN  INITIAL  OUTLET  TEMPERATURE. 

C  THE  MATRIX  RUNS  FROM  N=2  TO  N=NMAX-1. 

T01=0.0 

DO  37  J= (11-20) , (II-l) 

T0I=T0I+TG0 (J) 

37  CONTINUE 

T0I=T0I/20. 

TLI=0.0 

DO  36  J=(II-20) , (II-l) 

TLI=TLI+TGL(J) 

36  CONTINUE 

TLI=TLI/20. 

WRITE (53,*)  '  T0I=  ',T0I, '  TLI=  ',TLI 
DO  26  N=2,NMAX-1 

TRP  (N)  =  (TLI-TOI)  /  (NMAX-3)  *  (N-2)+T0I 
TGP  (N)  =TRP  (N) 

TTP  (N)=TRP  (N) 

26  CONTINUE 

C 


C  SET  THE  INITIAL  CONDITIONS  AT  THE  END  POINTS. 

TRP (1)=T0I 
TOP (1)=T0I 
TTP (1) =T0I 
TRP (NMAX)=TLI 
TOP (NMAX) =TLI 
TTP (NMAX)=TLI 
C 

C  DETERMINE  THE  REYNOLDS  NUMBER  AND  NUSSELT  NUMBER. 

REP S=RF*P S *DOTM/AF /ENU/RHO 
HNU=HTC*DH/EKG 

WRITE (53,*)  '  DOTM=  ',DOTM, '  REPS=  ',REPS,'  HNU=  '  ,  HNU 
WRITE(53,*)  DP,AF,RH0,  '  DP,AF,RHO' 

C 

C  DEFINE  SOME  OTHER  CONSTANTS  NEEDED  BELOW. 

TA=2 . 0*AS*HTC/DOTM/CP 
TB=2 . 0*RHO*AF*C/DOTM/GAM 
TC=HTC*AS/C/CM/RHOM/ (A* (1-POR) ) 

TJ=2 . 0*HTC*AST/DOTM/CP 

TD=TB-TA-TJ 

TE=EKS/DELX/RHOM/CM/C 

TF=HTC*AST/RHOM/CM/ACXT/C 

TH=1.0-2.0*TE-TF 

WRITE (53,*)  TA,  TB,TC,TD,TE,TF,TH,TJ,  '  TA, TB, TC, TD, TE, TF,  TH,  TJ 
C 

C  START  THE  INTERPOLATION  COUNTER,  ICOUNT,  AND  SET  THE  INITIAL 

C  INTERPOLATION  MARKER, 10. 

I0=II 

ICOUNT=l 

TG2 (I0)=TGP (NMAX) 

C 

27  DO  31  I=2,ITRAN 

C 

IF  (ICOUNT.EQ.IINT)  THEN 
10=10+1 
ICOUNT=0 
END  IF 
C 

C  INTERPOLATE  THE  VALUE  OF  THE  INCOMING  GAS  TEMPERATURE. 

TEMP=( (TGO (I0+1)-TG0 (10) )/IINT) *ICOUNT+TGO (10) 

C 

C  SET  THE  TEMPERATURE  OF  THE  GAS  AT  THE  INLET  TO  THE  MATRIX 

C  EQUAL  TO  THE  INTERPOLATED  EXPERIMENTAL  TGO. 

TGC(I)=TEMP 

C 

C  DETERMINE  THE  VALUE  OF  TG  AND  TM  AT  THE  NEXT  LOCATION. 

DO  22  N=2,NMAX-1 
C 

TGC  (N)  =  (TGP  (N-1)  -TGP  (N+1)  +TA*TRP  (N)  +TJ*TTP  (N)  +TD*TGP  (N) )  /TB 
TRC  (N)  =TC*TGP  (N)  +TRP  (N)  *  (1 . 0-TC) 

TTC  (N)  =TE*TTP  (N+1)  +TE*TTP  (N-1)  +TF*TGP  (N)  +TH*TTP  (N) 

C 

22  CONTINUE 

C 

C  SET  THE  TEMPERATURE  OF  THE  GAS/MATRIX  AT  THE  EXIT. 

TGC (NMAX) =TGC (NMAX-1) 

TRC  (NMAX)  =TRC  (NMAX-1 ) 

TTC (NMAX) =TTC (NMAX-1 ) 

C 

C  PUT  AN  ADIABATIC  CONDITION  AT  THE  FRONT  OF  THE  TUBE  AND  MATRIX 

TRC(1)=TRC(2) 

TTC(1)=TTC(2) 

C 

C  SET  THE  CURRENT  TIME  STEP  TEMPS  TO  THE  PREVIOUS  ONES. 

DO  55  J=1,NMAX 
TGP(J)=TGC(J) 


TRP  (J)=TRC(J) 

TTP (J)=TTC(J) 

55  CONTINUE 

C 

IF (ICOUNT.EQ.O)  THEN 
TG2 (I0)=TGC (NMAX) 

IF (ABS (TG2 (10) ) .GT.IOO.)  THEN 

WRITE (*,*)  'STABILITY  PROBLEM  OCCURING' 

STOP 
END  IF 

END  IF 
C 

I COUNT= I COUNT+ 1 
C 

31  CONTINUE 

C 

C 

C  CHECK  FOR  THE  MAXIMUM  SLOPE  AND  ELOC. 

SLOPE1=0.0 

DO  71  J=(II+5)  ,  (IF-5) 

SLOPE= (TG2 ( J-2) +8 . 0*TG2 ( J+l)-8 . 0*TG2 (J-1) -TG2 ( J+2) ) /12 . 0/SAMPT 
IF (SLOPE. GE. SLOPE!)  THEN 
SLOPEl=SLOPE 
LOC=J 
END  IF 

71  CONTINUE 

ELOC= (LOC-1) *SAMPT 

WRITE (53,*)  SLOPEl,  '  SLOPEl' ,  ELOC,  '  ELOC  ' 

C 

C  DETERMINE  THE  SPONGE  EFFECT  DELAY  TIME,  SEDT. 

TG2AV=TLI 
DO  83  J=II,IF 

TG2AV= (TG2AV* ( J-II+20) +TG2 (J) ) / (J-II+21) 

IF ( (TG2 (J+50)-TG2AV) .GE.0.4)  THEN 
SEDT=(J-II) *SAMPT 
TAU2=TAU0+SEDT 
GO  TO  84 
END  IF 

83  CONTINUE 

84  CONTINUE 
C 

WRITE (53,*)  '  SEDT  ',  SEDT,  ' ****TAU2  Equals',  TAU2 
C  OUTPUT  THE  TG  AND  TRs 

DO  57  J=II,IF 
TIME=(J-1) *SAMPT 
WRITE (52, 106)  TIME, TG2 ( J) , TGL ( J) 

57  CONTINUE 

C 

C  WRITE (53,  107)  MESH, REPS, HNU, RF 

C  OUTPUT  IMPORTANT  PARAMETERS 

WRITE  (53,*)  '  EL  ,  POR  ,  DELX  ,  MESH  ' 

WRITE (53,*)  EL, POR, DELX, MESH 
C 

101  FORMAT (F15 . 8, lx,F15 . 8, lx, F15 . 8, lx,  F15 . 8) 

102  FORMAT (F15.8/F15. 8) 

103  FORMAT (F15. 8) 

104  F0RMAT(F15.8,1X,F15.8) 

106  FORMAT (F15. 8, 1X,F15.8, lx,F15.8) 

107  FORMAT (4F15. 8) 

STOP 

END 
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PROGRAM  MAIN3C 

THIS  IS  A  COMPUTER  PROGRAM  TO  DETERMINE  THE  HEAT 
TRANSFER  COEFFICIENT  AND  FRICTION  FACTOR  FOR  A  WIRE-SCREEN 
REGENERATOR  SUBJECTED  TO  A  QUASI-STEP  CHANGE 
IN  INLET  TEMPRATURE. 

THE  WORKING  FLUID  IS  HELIUM  AT  LARGE  PRESSURES  (5-20  ATM) . 

THE  MASSFLOW  RATE  IS  GIVEN  (MEASURED)  AND  ASSUMED  STEADY. 

A  FINITE  DIFFERENCE  SOLUTION  TO  THE  ENERGY  EQUATION  FOR 
INCOMPRESSIBLE  FLOW  IS  USED  TO  CALCULATE  THE  TEMPERATURE 
OF  THE  GAS  AND  MATRIX  AS  A  FUNCTION  OF  TIME  AND  SPACE. 

THE  VALUES  OF  THE  HTC  ARE  DETERMINED  ITERATIVELY  BY  MATCHING 
THREE  CRITERIA; 

1)  THE  SIZE  AND  LOCATION  OF  THE  MAXIMUM  SLOPE  OF  THE 
TEMP-TIME  CURVE  OF  THE  TRANSIENT. 

2)  THE  TIME  IT  TAKES  FOR  THE  OUTLET  GAS  TEMPERATRUE 

TO  DEVIATE  SIGNIFICANTLY  FROM  THE  INITIAL  REGENERATOR 
TEMPERATURE,  THE  SPONGE  EFFECT  DELTA  T,  SEDT . 

3)  THE  ROOT-MEAN-SQUARE  DIFFERENCE  OF  THE  EXPERIMENTAL 
AND  CALCULATED  OUTLET  TEMPERATURES. 

A  FUNCTION  SUBPROGRAM  NAMED  ZEROIN  IS  USED  TO  FIND  THE 
SOLUTION  OF  F(X)=0.  ANOTHER  FUNCTION  SUBPROGRAM, 

FUNSMS/DT/RM,  IS  USED  TO  DETERMINE  THE  MAXIMUM  SLOPE,  SEDT,  AND 
ROOT-MEAN-SQUARE-DIFFERENCE,  RMSD,  OF  THE  CALCULATED 
OUTLET  TEMPERATURE  PROFILE,  TG2,  FROM  THE  EXPERIMENTAL  ONE, 
GIVEN  THE  INPUT  TEMPERATURE-TIME  TRACE  AND  A  GUESS  AT 
THE  HEAT  TRANSFER  COEFFICIENT,  HTC. 

FUNSMS/DT/RM  IS  A  NUMERICAL  SOLUTION  FOR  THE  TRANSIENT 
TEMPERATURE  DISTRIBUTION  FOR  STEADY  INCOMPRESSIBLE 
FLOW  THROUGH  A  WIRE-SCREEN  REGENERATOR  SUBJECTED  TO 
A  QUASI-STEP  CHANGE  IN  TEMPERATURE  BETWEEN  TOI  AND  TLI . 

DOTM  IS  THE  MASS  FLOW  RATE  (FUNCTION  OF  TIME  ONLY)  . 

TGC(N)  ,TRC(N)  ,AND  TTC  (N)  ARE  THE  GAS,  REGENERATOR,  AND  TUBE 
TEMPERATURES  DURING  THE  CURRENT  TIME-STEP  OF  THE  TRANSIENT. 

TGP  (N)  ,TRC(N)  ,  AND  TTC(N)  ARE  THE  GAS,  REGENERATOR,  AND  TUBE 
TEMPERATURES  DURING  THE  PREVIOUS  TIME-STEP. 

TGO(I)  AND  TGL(I)  ARE  THE  EXPERIMENTALLY  MEASURED  GAS 
TEMPERATURES  AT  THE  INLET  AND  EXIT  RESPECTIVELY. 

TG2(I)  IS  THE  CALCULATED  OUTLET  TEMPERATURE. 

TAV,  PAV  ARE  THE  EXPERIMENTALLY  MEASURED  AVERAGE  TEMPERATURE 
AND  PRESSURE  OF  THE  GAS  DURING  THE  TRANSIENT. 

D  IS  THE  PIPE  DIAMETER;  EL  IS  THE  REGENERATOR  LENGTH. 

POR  IS  THE  PORISITY;  I,J,N  ARE  THE  TIME  AND  SPACE  COUNTERS. 

DH  IS  THE  HYDRAULIC  DIAMETER  OF  THE  REGENERATOR  MESH. 

DELT,  DELX  ARE  THE  TIME  AND  SPACE  INCREMENTS  DURING  THE  ACTUAL 
TRANSIENT. 

JMAX  AND  NMAX  ARE  THE  MAXIMUM  TIME  AND  SPACE  INCREMENTS  OF 
THE  TOTAL  DATA  SAMPLE. 

ITRAN  IS  THE  NUMBER  OF  TIME  INCREMENTS  DURING  THE  ACTUAL 
TRANSIENT,  (TAUF-TAUO) /DELT . 

10  IS  THE  COUNTER  FOR  THE  EXPERIMENTAL  DATA  AT  THE  BEGINING 
OF  THE  CURRENT  I COUNT. 

IINT  IS  THE  NUMBER  OF  CALCULATION  TIME  INCREMENTS  BETWEEN 
EXPERIMENTAL  DATA  POINTS. 

ICOUNT  IS  THE  COUNTER  FOR  THE  CALCULATION  SUBINTERVALS 
BETWEEN  EXPERIMENTAL  DATA  POINTS. 

11  AND  IF  ARE  THE  TIME  INDEX  OF  THE  START  AND  END  OF  THE  ACTUAL 
TRANSIENT  ASSOCIATED  WITH  II  AND  IF. 

EMU  IS  THE  DYNAMIC  VISCOUSITY;  ENU  THE  KINEMATIC  VISCOUSITY. 

REPS  IS  THE  REYNOLDS  NUMBER  BASED  ON  PORE  SIZE. 

REDW  IS  THE  REYNOLDS  NUMBER  BASED  ON  WIRE  DIAMETER. 

CD  IS  A  DRAG  COEFFICIENT  USED  IN  F3C . 

PS  IS  THE  DISTANCE  BETWEEN  WIRES  IN  THE  SCREENS. 

PSC  IS  THE  ACTUAL  PORE  SIZE;  PSC=  (PS*PS*2*DP*RF)  **0 . 33333  . 

DP  IS  THE  WIRE  DIAMETER. 

PIT  IS  THE  PITCH  =  DP+PS. 

HTC  IS  DEFINED  AS  Q  /  A/  TG-TR 


C  R  IS  THE  GAS  CONSTANT  FOR  THE  WORKING  FLUID. 

C  EKG  IS  THE  THERMAL  CONDUCTIVITY  OF  THE  GAS. 

C  EKS  IS  THE  THERMAL  CONDUCTIVITY  OF  THE  304  STAINLESS  STEEL. 

C  A  IS  THE  CROSS  SECTIONAL  AREA  OF  THE  PIPE. 

C  AF  IS  THE  FLOW  AREA  OF  THE  REGENERATOR;  AF  =  POR*A. 

C  ASP  IS  THE  SURFACE  AREA  FOR  HEAT  TRANSFER  PER  UNIT  VOLUME. 

C  AS  IS  THE  HEAT  TRANSFER  AREA  PER  SPACE  INCREMENT. 

C  ATOT  IS  THE  TOTAL  HEAT  TRANSFER  AREA  IN  THE  SCREEN  UNIT. 

C  AST  IS  THE  SURFACE  AREA  OF  THE  TUBE  FOR  HEAT  TRANSFER. 

C  ACXT  IS  THE  CROSS  SECTIONAL  AREA  OF  TUBE  WALL  FOR  CONDUCTION. 

C  CP  IS  THE  SPECIFIC  HEAT  OF  THE  GAS  AT  CONSTANT  PRESS. 

C  CV  IS  THE  SPECIFIC  HEAT  OF  THE  GAS  AT  CONSTANT  VOLUME. 

C  GAM  IS  THE  RATIO  OF  CP/CV. 

C  HTCMS  IS  THE  HEAT  TRANSFER  COEFFICIENT  DETERMINED  BY 

C  USING  THE  MAXIMUM  SLOPE  CRITERIA. 

C  HTCDT  IS  THE  HEAT  TRANSFER  COEFFICIENT  DETERMINED  BY 

C  USING  THE  SPONGE  EFFECT  DELAY  TIME  CRITERIA. 

C  HTCRM  IS  THE  HEAT  TRANSFER  COEFFICIENT  DETERMINED  BY 

C  USING  THE  ROOT  MEAN  SQUARE  DIFFERENCE  CRITERIA. 

C  HTCMAX, HTCMIN  ARE  BOUNDING  GUESSES  FOR  HTC  IN  ZEROIN. 

C  HTCSTAB  IS  A  BOUNDING  HTC  FOR  STABILITY. 

C  TOL  IS  A  TOLERANCE  USED  IN  ZEROIN. 

C  ESLOPE  AND  ELOC  ARE  THE  MAGNITUDE  AND  LOCATION  IN  TIME  OF 

C  THE  MAXIMUM  SLOPE  OF  THE  EXPERIMENTAL  TIME-TEMP  CURVE. 

C  LOC  IS  THE  INDEX  OF  ELOC. 

C  CSLOPE  IS  THE  MAXIMUM  SLOPE  OF  THE  TEMPERATURE-TIME  CURVE 

C  AT  TIME  EQUALS  ELOC  CALCULATED  FROM  THE  MODEL. 

C  RHO  IS  THE  GAS  DENSITY. 

C  CM  IS  THE  MATRIX  HEAT  CAPACITY;  RHOM  THE  MATRIX  DENSITY. 

C  THK  IS  THE  TUBE  WALL  THICKNESS;  WSU  THE  MASS  OF  THE  SCREEN  UNIT. 

C  SCTHK  IS  THE  THICKNESS  OF  A  SCREEN. 

C  MESH  IS  THE  MESH  OF  THE  SCREENS  (PER  INCH;  SEE  MIYABE.) 

C  SAMPT  IS  THE  EXPERIMENTAL  DATA  SAMPLING  TIME. 

C  TIME  IS  THE  TIME  DURING  THE  CALCULATED  SUBINTERVAL. 

C  TTOT  IS  THE  TOTAL  DATA  SAMPLING  TIME,  I.E.  JMAX=TTOT/SAMPT . 

C  TAUO  AND  TAUF  ARE  THE  BEGINNING  AND  ENDING  TIMES  FOR  THE 

C  ACTUAL  TRANSIENT. 

C  TAU2  IS  THE  TIME  AT  THE  END  OF  THE  SEDT;  12  IS  ITS  INDEX. 

C  TO  I  AND  TLI  ARE  THE  INITIAL  TEMPERATURES  AT  THE  INLET  AND  OUTLET. 

C  THETA  IS  THE  ANGLE  OF  CONTACT  BETWEEN  STRANDS  OF  THE  SCREEN. 

C  RF  IS  THE  REDUCTION  FACTOR  FOR  THE  SCREEN;  RF=EL  (ACT) /ELO  . 

C  RMSD  IS  A  ROOT  MEAN  SQUARE  DIFFERENCE. 

C  DPSU  IS  THE  AVERAGE  EXPERIMENTALLY  MEASURED  PRESSURE  DROP 

C  THROUGH  THE  REGENERATOR. 

C  F1/2/3E  AND  F1/2/3C  ARE  EXPERIMENTAL  AND  CALCULATED  FRICTION 

C  FACTORS . 

C  USTAR  IS  A  VELOCITY  BASED  ON  THE  RATIO  OF  FRONTAL  AREA  FOR 

C  FLOW  TO  THE  TOTAL  AREA. 

C  BETA  IS  THE  RATIO  OF  THE  MINIMUM  FLOW  AREA  TO  THE  TOTAL  AREA. 

C  ARP  IS  AN  AREA  FACTOR  USED  FOR  FRICTION  FACTOR  F2CARP. 

C  ELARM  IS  A  LENGTH  DEFINED  BY  ARMOUR. 

C  REARM  IS  A  REYNOLDS  NUMBER  BASED  ON  ELARM. 

C  F2CARP  IS  ARMOUR'S  FRICTION  FACTOR  FOR  ROUND  WIRE  SCREENS. 

C 

DOUBLE  PRECISION  TGO  (20001) , TGL  (20001) , TG2  (20001) 

INTEGER  I,N,NMAX, UMAX, LOC, II, IF, 12, IINT,  ITRAN,  J 
DOUBLE  PRECISION  R,  RHO,  CM,  RHOM,  A,  CP,  CV,  GAM,  DOTM,  ACXT,  ARP, 

$DELT,  DELX,  DH,  AF, ENU,  EMU, POR, EL,  D,  PI, DP, PS, EKG,  SAMPT,  SEDT, 

$TAV,  TOL,  AS,  C,  P IT,  TOI ,  TLI,  PAV,  HTC,  THK,  EKS,  WSU,  THETA,  HTCSTAB,  ASP, 
SESLOPE,  ELOC,  HTCMAX,  HTCMIN,  REPS,  TAUO,  TAUF,  MESH,  RF,  AST, 

$HTCMS,  HTCDT,  HTCRM,  SZERMS,  SZERDT,  SZERRM,  TAU2,FUNSMS,  FUNSDT,  FUNSRM, 
$DPSU,  FIE,  F2E,  F3E,  FlC,  F2C,  F3C,  USTAR,  Q,  ATOT,  SCTHK,  BETA,  CD,  REDW, 
$REARM, ELARM, F2CARP , PSC 
C 

COMMON  TGO,  TGL,  DOTM,  AS,  CP,  DELX,  SAMPT, RHO,  AF,  C,  CM,  GAM, DELT,  THK, 
$POR,  RHOM,  TOI,  TLI,  DP,  ENU,  EKG,  DH,  PS,  II,  IF,  IINT,  ITRAN, 


oo  oo  o  nooo  oo  oo  oo  ooo  oo  ooo  o  non  on  ooo  ono  ooo  oo 


$D, PI, A, ACXT, AST,  EKS, NMAX, LOG 
C 

EXTERNAL  SZERMS,  SZERDT,  SZERRM,  FLTNSMS,  FUNSDT,  FUNSRM 
INTRINSIC  COS, MAX, SQRT, ABS,NINT, ATAN 

OPEN  THE  DATA  INPUT/OUTPUT  FILES. 

OPEN(50,  FILE=' mss.  in',  STATUS=' UNKNOWN') 

OPEN (52,  FILE=' mss. out', STATUS=' UNKNOWN') 

OPEN(53,  FILE=' final.dat' ,  STATUS=' UNKNOWN' ) 

INPUT  FLOW  PARAMETERS  AND  NUMERICAL  VARIABLES. 

INPUT  THE  WIRE  DIAMETER  AND  PORE  SIZE. 

READ (50, 102)  DP, PS 
WRITE (*,*)  'DP  =' ,  DP,  'PS  =  '  ,  PS 

INPUT  REGENERATOR  LENGTH  AND  WEIGHT. 

READ (50, 102)  EL,WSU 

WRITE  (*,*)  'LENGTH  =  ',EL, '  WEIGHT  =  ',  WSU 

INPUT  THE  DELT  FOR  THE  TRANSIENT. 

WRITE (*,*)  '  INPUT  DELT  ' 

READ(*,*)  DELT 

INPUT  THE  MESH  NUMBER 
READ (50, 103)  MESH 
WRITE (*,*)  'MESH  =  ',  MESH 

THE  TUBE  DIAMETER  AND  THICKNESS  ARE  CONSTANT. 

D=0. 014859 
THK=0. 000508 

WRITE  (*,*)  'INPUT  THE  MAXIMUM  SPACE  INCREMENTS' 

READ(*,*)  NMAX 
WRITE  (*,*)  'NMAX  =  ',  NMAX 

R  IS  THE  GAS  CONSTANT  FOR  HELIUM  IN  J/Kg*K 
R=2077. 

DETERMINE  THE  SPACE  INCREMENTS. 

DELX=EL/ (NMAX-1) 

WRITE (*, *)DELX,  '  DELX' 

ADD  TWO  SPACE  INCREMENTS  TO  THE  TOTAL. 

NMAX=NMAX+2 

EKG  IS  THE  THERMAL  CONDUCTIVITY  OF  HELIUM.  (SEE  KAYS)  . 
EKG=0.155 

EKS  IS  THE  THERMAL  CONDUCTIVITY  OF  304  SS,  W/m-K.  (INCROPERA) 
EKS=14.9 

NEED  A  VALUE  FOR  EMU  FOR  THE  HELIUM  (Pa*SEC)  . 

EMU=1 . 95E-05 


RHO  IS  THE  DENSITY  OF  THE  GAS,  ASSUMED  CONSTANT,  BASED  ON  THE 
AVERAGE  TMEPERATURE  AND  PRESSURE  AT  THE  VENTURI  DURING  THE 
TRANSIENT. 

READ  (50,  102)  PAV,  TAV 

WRITE (*,*)  'PAV  =  ',  PAV,  '  TAV  =  ',  TAV 
RHO= (PAV* 6895 .) /R/ (TAV+273 . ) 

DETERMINE  THE  KINEMATIC  VISCOSITY. 

ENU=EMU/RHO 

DEFINE  PI 
PI=4.0*ATAN(1. ) 


oooo  ooooooooo  M  ooooo  ooo  nonoooooooooo 


NEED  THE  CROSS-SECTIONAL  AREA  OF  THE  PIPE. 

A=PI*D**2/4.0 

NEED  THE  CROSS-SECTIONAL  AREA  OF  THE  TUBE  WALL. 

ACXT=PI* (D+THK/2 .0) *THK 

NEED  A  VALUE  FOR  THE  THERMAL  CAPACITY  OF  THE  MATRIX 
MATERIAL  IN  J/KG/K. 

CM=477.0 

NEED  A  VALUE  FOR  THE  DENSITY  OF  THE  MATRIX  IN  KG/M3 . 

RHOM=7100 . 

NEED  A  VALUE  FOR  THE  HEAT  CAPACITY  OF  THE  GAS  IN  J/KG/K. 

CP=5193. 

NEED  A  VALUE  FOR  GAM  AND  CV. 

CV=CP-R 

GAM=CP/CV 

NEED  A  VALUE  FOR  HEAT  TRANSFER  AREA,  AS.  (SEE  MY  NOTES) 

ENTER  THE  REDUCTION  FACTOR. 

READ  (50,  103)  RF 
PIT=DP+PS 

DETERMINE  THE  CONTACT  ANGLE,  THETA. 

THETA=ATAN (DP/PIT) 

WRITE (*,*)  '  THETA,  PIT,  RF  ',  THETA,PIT,RF 

DETERMINE  THE  SURFACE  AREA  PER  UNIT  ORIGINAL  VOLUME. 

ASP=PI* (SQRT(PIT**2+DP**2)-(DP*THETA**2/PI) )  /PIT**2 
CORRECT  THE  SPECIFIC  AREA  FOR  REDUCTION  FACTOR.  (SEE  NOTES) 

IF  (RF.LT.0.65)  THEN 
ASP=ASP*0.8 
GO  TO  15 

ELSEIF  (RF.LT.0.8)  THEN 
ASP=ASP* 0.970 
GO  TO  15 

ELSEIF  (RF.LT.0.95)  THEN 
ASP=ASP* 0.9992 
GO  TO  15 

ELSE 

GO  TO  15 

END  IF 
CONTINUE 

AS=ASP*A* (DELX/RF) 

ATOT=AS*EL/DELX 

WRITE (*,*)  '  THE  TOTAL  HEAT  TRANSFER  AREA  IS  ',  ATOT 

DETERMINE  THE  POROSITY  (SEE  MIYABE  1982) . 

POR=l . 0- (WSU/RHOM/A/EL) 

WRITE (*,*)  POR,  '  THE  POROSITY  IS' 

DETERMINE  THE  FLOW  AREA. 

AF=POR*P I *D * *2 / 4 . 0 

NEED  THE  TUBE  HEAT  TRANSFER  AREA. 

AST=PI*D*DELX*POR 

DEFINE  THE  HYDRAULIC  DIAMETER. 

DH=EL*AF/ATOT 

WRITE (*,*)  '  DH  =',DH 

FIND  THE  NUMBER  OF  DATA  POINTS  PER  TRANSIENT,  UMAX. 

INPUT  THE  SAMPLING  TIME  AND  TOTAL  NUMBER  OF  EXPERIMENTAL  POINTS. 


READ (50, 102)  SAMPT, TTOT 
JMAX=INT (TTOT /SAMPT) 

C  WRITE (*,*)  JMAX,  '  JMAX' 

C 

C  DEFINE  THE  COURANT  NUMBER. 

C=DELX/DELT 

C  WRITE  (*,*)  C,  '  THE  COURANT  NUMBER  IS  ' 

C 

N=1 

1=1 

J=1 

C 

C  ENTER  THE  EXPERIMENTALLY  MEASURED  VALUES  OF  DOTM  AND  THE 

C  GAS  TEMPERATURES  WHERE  TGO  APPROXIMATES  A  STEP  CHANGE. 

READ (50, 103)  DOTM 

DOTM  =  l.l*DOTM 

WRITE(*,*)  '  DOTM  =  ',DOTM 

DO  21  J=1,JMAX 

READ (50, 104)  TGO ( J) , TGL ( J) 

21  CONTINUE 

C 

C  DETERMINE  THE  INDEX  OF  THE  BEGINNING  OF  THE  ACTUAL  TRANSIENT 

C  II,  AND  THE  END  OF  THE  ACTUAL  TRANSIENT,  IF . 

C 

READ (50, 102)  TAU0,TAUF 
II=NINT (TAUO /SAMPT) +1 
IF=NINT (TAUF/SAMPT) +1 

C  WRITE (*,*)  TAUO, TAUF, II, IF, ' ***********  TAUO, TAUF, II, IF  ' 

C 

C  DETERMINE  THE  INDEX  FOR  THE  END  OF  THE  SEDT. 

READ (50, 103)  TAU2 
I2=NINT (TAU2/SAMPT) +1 

C  WRITE (*,*)  'TAU2  AND  12  ARE  ',  TAU2, 12 

SEDT=(I2-II) *SAMPT 
C 

C  DETERMINE  THE  NUMBER  OF  INCREMENTS  BETWEEN  DATA  POINTS,  IINT,  AND 

C  THE  NUMBER  OF  TIME  INCREMENTS  DURING  THE  ACTUAL  TRANSIENT,  ITRAN. 

IINT=NINT (SAMPT/DELT) 

ITRAN=NINT ( (TAUF-TAUO) /DELT) +1 
C  WRITE (*,*)  IINT, ITRAN, '  IINT, ITRAN' 

C 

C  DETERMINE  THE  INITIAL  MEAN  TEMPERATURE  AT  BOTH  THE  OUTLET 

C  AND  INLET,  TLI  AND  TOI. 

TLI=0.0 

DO  37  J= (11-20) , (II-l) 

TLI=TLI+TGL (J) 

37  CONTINUE 

TLI=TLI/20. 

C 

T0I=0.0 

DO  36  J= (11-20), (II-l) 

T0I=T0I+TG0 (J) 

36  CONTINUE 

T0I=T0I/20. 

WRITE (*,*)  TLI, TOI,'  TLI, TOI  ' 

C 

C  ENTER  THE  MAXIMUM  EXPERIMENTALLY  MEASURED  SLOPE,  ESLOPE, 

C  OF  THE  OUTLET  GAS  AND  ITS  LOCATION  ON  THE  TIME-TEMP  CURVE. 

READ (50, 102)  ESLOPE, ELOC 
LOC=NINT (ELOC/SAMPT) +1 

C  WRITE(*,*)  ESLOPE,ELOC, LOC, '  ******ESLOPE, ELOC,  LOC' 

C 

C  ENTER  THE  EXPERIMENTALLY  MEASURED  VALUES  OF  THE  PRESSURE 

C  DROP  THROUGH  THE  REGENERATOR,  DPSU,  AND  THE  SCREEN  THICKNESS. 

READ (50, 102)  DPSU, SCTHK 

C  WRITE  (*^*)  '  THE  DPSU  IS  ',DPSU, '  SCTHK  IS  ',  SCTHK 


non  on  on  o  nononn  n  n  n  n  oonoo 


C 

c 

C  USE  ZERO  IN  TO  DETERMINE  THE  HTC  GIVEN  BOUNDING  VALUES  OF 

C  THE  HTC  AND  THE  LOCATION  AND  MAGNITUDE  OF  THE  EXPERIMENTALLY 

C  MEASURED  MAXIMUM  SLOPE  AND  SEDT  OF  THE  TRANSIENT. 

TOL=5 . OE+1 

SIGNIFY  THE  MAXIMUM  HTC  FOR  THIS  CONDITION. 
HTCSTAB=C*CP*RHO*AF/ (AS+AST) /GAM 

WRITE (*,*)  'THE  MAXIMUM  HTC  FOR  STABILITY  IS  ABOUT  ',HTCSTAB 

INPUT  THE  GUESSES  FOR  HTCMIN  AND  HTCMAX 
WRITE (*,*)  '  INPUT  HTCMIN  AND  HTCMAX  ' 

READ ( * , * )  HTCMIN, HTCMAX 

WRITE ( * , * )  '  HTCMIN=' , HTCMIN, '  HTCMAX=' ,  HTCMAX 

DETERMINE  THE  HTCs  BASED  ON  THE  THREE  CRITERIA. 

HTCRM=SZERRM  (HTCMIN,  HTCMAX,  FUNSRM,  TOL,  12 ) 

HTCMS=SZERMS  (HTCMIN,  HTCMAX,  FUNSMS,  TOL,  ESLOPE) 

HTCDT=SZERDT  (HTCMIN,  HTCMAX, FUNSDT,  TOL,  12,  TAU2) 

WRITE (*,*)  HTCDT, '  =HTCDT' 

DETERMINE  THE  FRICTION  FACTORS. 

CHANGE  THE  PRESSURE  DROP  TO  PASCALS. 

DPSU=DPSU*6895.0 

FI  IS  ACCORDING  TO  CHEN;  200  IS  THE  NUMBER  OF  SCREENS. 

BETA=(  1.0-DP /PIT) **2 
USTAR=DOTM/RHO/A/BETA 
WRITE (*,*)  '  USTAR  =  ',USTAR 
REPS=PS*USTAR/ENU 
F1E=2 . 0*DPSU/200 . /RHO/USTAR**2 
F1C=49.78/REPS**0. 968+0. 318 

F2  IS  DUE  TO  ARMOUR  AND  CANNON;  Q  IS  A  TORTUOUSITY  FACTOR. 

Q=l. 

PSC= (PS**2*2 . *DP*RF) ** (1 . /3  . ) 

ELARM=1 . 0/ (ASP**2*PSC) 

REARM=DOTM*ELARM/A/EMU 
ARP=PI*SQRT(PIT**2+DP**2) /PIT**2 
F2E=DPSU*POR*  *2  *PSC*A*  *2  *RHO/Q/EL/DOTM**2 
F2CARP=8 . 61*EMU*ARP**2*PSC*A/DOTM+0 . 52 
F2C=8 . 61/REARM+O . 52 

WRITE (*,*)  '  REARM  =  REARM, '  ELARM  =  ',ELARM 

F3  IS  THE  SAME  AS  WIESE  DUE  TO  KAYS/LONDON. 

REDW=DP  *USTAR/ENU 

CD=10.0**  (  (1.1/  (POR**2)/REDW**  (0 . 254) ) - (0 . 54/POR)  ) 

F3E=DPSU* (AF**3) *RHO*DP*RHOM/2 . 0/DOTM**2/WSU 
F3C=CD* (DH/SCTHK) * (POR/BETA) **2 
WRITE (*,*)  '  FIE  F2E 

WRITE (*,106)  F1E,F2E,F3E 
WRITE (*,*)  '  FlC  F2C 

WRITE(*,106)  F1C,F2C,F3C 
WRITE(*,*)  '  HTCDT  =  HTCDT 

OUTPUT  THE  RESULTS  TO  FILE. 

REP  S=P  S*DOTM/ AF /EMU 

WRITE(53, 107) 

WRITE  (53,  105)  HTCMS, ESLOPE, HTCDT,  SEDT 
WRITE (53, 108) 

WRITE (53, 105)  HTCRM, MESH, REPS, RF 

WRITE (53,*)  '  FIE  F2E  F3E' 

WRITE (53, 106)  F1E,F2E,F3E 


F3E  ' 
F3C' 


WRITE (53,*)  '  FlC  F2C  F3C' 

WRITE (53, 106)  F1C,F2C,F3C 
WRITE (53,*)  'DELT=  ',DELT, '  NMAX=  ',NMAX 
WRITE (53,*)  '  FOR  =  ',POR, '  DH=  ',DH 
C 

102  FORMAT (F15.8/F15. 8) 

103  FORMAT (F15. 8) 

104  F0RMAT(F15.8,1X,F15.8) 

105  FORMAT  (4F15. 8) 

106  FORMAT (F15. 8, 1X,F15.8, 1X,F15.8) 

107  FORMAT  (5X,  '  HTCMS'  ,  9X,  '  ESLOPE' ,  lOX, '  HTCDT'  ,  lOX,  '  SEDT' ) 

108  FORMAT  (5X,  'HTCRM'  ,  IIX,  'MESH' ,  IIX,  'REPS' ,  12X,  'RF'  ) 

STOP 

END 

C 

C 

C  THIS  SUBROUTINE  FINDS  THE  VALUE  OF  HTC  WHICH  BEST  FITS  THE 

C  VALUE  OF  ESLOPE. 

C 

DOUBLE  PRECISION  FUNCTION  SZERMS  (HTCMIN,  HTCMAX,  FUNSMS,  TOL, 
$ESLOPE) 

DOUBLE  PRECISION  HTCMIN, HTCMAX, FUNSMS, TOL,  ESLOPE 
C 

C  A  ZERO  OF  THE  FUNCTION  F  (X)  IS  COMPUTED  IN  THE  INTERVAL  AX,BX 

C 

C  INPUT 
C 

C  AX  LEFT  ENDPOINT  OF  INITIAL  INTERVAL 
C  BX  RIGHT  ENDPOINT  OF  INITIAL  INTERVAL 

C  F  FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F(X)  FOR  ANY  X  IN 

C  THE  INTERVAL  AX,BX 

C  TOL  DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THE 
C  FINAL  RESULT 

C 
C 

C  OUTPUT 
C 

C  ZEROIN  ABCISSA  APPROXIMATING  A  ZERO  OF  F  IN  THE  INTERVAL  AX,BX 
C 

C  IT  IS  ASSUMED  THAT  F  (AX)  AND  F  (BX)  HAVE  OPPOSITES  SIGNS 

C  WITHOUT  A  CHECK.  ZEROIN  RETURNS  A  ZERO  X  IN  THE  GIVEN  INTERVAL 

C  AX,BX  TO  WITHIN  A  TOLERANCE  4*MACHEPS*ABS (X)  +  TOL,  WHERE  MACHEPS 
C  IS  THE  RELATIVE  MACHINE  PRECISION. 

C  THIS  FUNCTION  SUBPROGRAM  IS  A  SLIGHTLY  MODIFIED  TRANSLATION  OF 
C  THE  ALGOL  60  PROCEDURE  ZERO  GIVEN  IN  R.  BRENT,  ALGORITHMS  FOR 

C  MINIMIZATION  WITHOUT  DERIVATIVES,  PRENTICE -HALL,  INC.  (1973) 

C 

C 

DOUBLE  PRECISION  A,  B, C, D,E, EPS, FA, FB, FC, TOLl,  XM,  P, Q,R, S 
C 

EXTERNAL  FUNSMS 
C 

C  COMPUTE  EPS,  THE  RELATIVE  MACHINE  PRECISION 
C 

EPS  =1.0 

61  EPS=EPS/2.0 

TOLl  =  1.0  +  EPS 
IF (TOLl. GT. 1.0) GO  TO  61 
C 

C  INITIALIZATION 
C 

A=HTCMIN 

B=HTCMAX 

FA=FUNSMS  (A,  ESLOPE) 

FB=FUNSMS  (B, ESLOPE) 


WRITE (*,*)  'fa=  ',FA, '  fb=  ',FB 
C 

IF( (FA/FB) .GT.0.0)  GO  TO  71 
C 

C  BEGIN  STEP 

C 

62  C=A 
FC=FA 
D=B-A 
E=D 

63  IF  (ABS (FC)  -GE.ABS (FB) )GO  TO  64 

A=B 

B=C 

C=A 

FA=FB 

FB=FC 

FC=FA 

C 

C  CONVERGENCE  TEST 
C 

64  TOLl  =  2.0*EPS*ABS (B)  +  0.5*TOL 

XM  =  0.5* (C-B) 

IF(ABS(XM)  .LE.TOLDGO  TO  69 
IF(FB.EQ.O.O)GO  TO  69 
C 

C  IS  BISECTION  NECESSARY? 

C 

IF  (ABS (E)  .LT. TOLl) GO  TO  67 
IF  (ABS (FA)  .LE.ABS (FB) )GO  TO  67 
C 

C  IS  QUADRATIC  INTERPOLATION  POSSIBLE? 

C 

IF(A.NE.C)GO  TO  65 
C 

C  LINEAR  INTERPOLATION 
C 

S=FB/FA 
P=2.0*XM*S 
Q=1.0-S 
GO  TO  60 
C 

C  INVERSE  QUADRATIC  INTERPOLATION 
C 

65  Q=FA/FC 
R=FB/FC 
S=FB/FA 

P=S* (2 . 0*XM*Q* (Q-R) - (B-A) * (R-1 . 0) ) 

Q  =  (Q  -  1.0)*(R  -  1.0)*(S  -  1.0) 

C 

C  ADJUST  SIGNS 
C 

60  IF(P.GT.0.0)Q  =  -Q 

P  =  ABS(P) 

C 

C  IS  INTERPOLATION  ACCEPTABLE? 

C 

IF(  (2.0*P)  .GE.  (3.0*XM*Q  -  ABS  (TOLl*Q) )  )GO  TO  67 
IF(P.GE.ABS(0.5*E*Q) )GO  TO  67 
E=D 
D=P/Q 
GO  TO  68 
C 

C  BISECTION 
C 

67  D=XM 

E=D 


c 

C  COMPLETE  STEP 
C 

68  A=B 
FA=FB 

IF(ABS(D) .GT.TOLl)  B  =  B  +  D 

IF (ABS (D) .LE.TOLl)  B  =  B  +  SIGN (TOLl, XM) 

FB  =  FUNSMS (B,E SLOPE) 

WRITE (*,*)  '  New  Value  of  B  ',B 
IF ( (FB* (FC/ABS (FC) ) ) .GT.O.O)GO  TO  62 
GO  TO  63 
C 

C  DONE 
C 

69  SZERMS  =  B 
RETURN 

71  WRITE (*,*)  ' FUNSMS  (BAX) /FUNSMS (BBX)  .GT. 0.0' 

STOP 

END 

C 

Q  ★★*★•*■*•*••*'•*■**★•*■*★**★*■*•****■*•★•*•****•*****★★*★★**★*********★★**★** 

C 

C  THIS  FUNCTIOM  SUBPROGRAM  DETERMINES  THE  VALUE  OF 

C  ESLOPE-CSLOPE  FOR  A  GIVEN  HTC . 

C 

DOUBLE  PRECISION  FUNCTION  FUNSMS (HTC, ESLOPE) 

C  THIS  IS  A  FUNCTION  SUBROUTINE  WHICH  DETERMINES  THE  VALUE 

C  OF  THE  DIFFERENCE  BETWEEN  THE  EXPERIMENTAL  AND  CALCULATED 

C  MAXIMUM  SLOPE  OF  THE  TIME-TEMPERATURE  CURVE  TRANSIENT. 

C  IT  IS  USED  TO  GET  THE  BEST  VALUE  FOR  HTC  IN  CONJUNCTION 

C  WITH  THE  FUNCTION  SZERMS (HTCMIN, HTCMAX, FA/FB,  TOL)  . 

C  msie.f  IS  THE  CONTROLING  PROGRAM  UNIT.  SEE  PROGRAM  tsie.f 

C  FOR  OTHER  DETAILS. 

C 

C  TA,  TB, TC,TD,TE,TF, TH,TI  ARE  CONSTANTS  USED  IN  CALCULATIONS. 

C  CSLOPE  IS  THE  CALCULATED  MAXIMUM  SLOPE. 

C  ESLOPE  IS  THE  EXPERIMENTALLY  MEASURED  SLOPE. 

C  LOC  IS  THE  LOCATION  OF  THE  MAX  SLOPE. 

C 

C 

DOUBLE  PRECISION  TGC (155) , TRC (155) , TTC (155) , TGP (155) , TRP (155) , 
$TTP (155) ,  TGO (20001) , TGL (20001) , TG2 (20001) 

INTEGER  I,N,NMAX, II, IF, LOC, 10, ITRAN, IINT, ICOUNT,  J 
DOUBLE  PRECISION  RHO, TA, TB, TC, TD, CM, RHOM, A, CP,  GAM, 

$DELT,  DELX, DH, AF, ENU, POR, D, PI, DP, PS, EKG, EKS,  HNU, 

$TE,  TF,  TH,  TJ, TOI, TLI, CSLOPE, ESLOPE, SAMPT, 

$AS , C, HTC, DOTM, REPS , TEMP , THK, ACXT, AST 
C 

COMMON  TGO,  TGL,  DOTM,  AS,  CP,  DELX,  SAMPT,  RHO,  AF,  C,  CM,  GAM,  DELT,  THK, 
$POR,  RHOM,  TOI,  TLI,  DP,  ENU,  EKG,  DH,  PS,  II,  IF,  IINT,  ITRAN, 
$D,PI,A,ACXT,AST,EKS,NMAX,LOC 
C 

INTRINSIC  COS, MAX, SORT, ABS, NINT,ATAN 
C  AT  THE  INITIAL  TIME,  THE  REGENERATOR  AND  GAS  ARE  AT  A 

C  STEADY  STATE  INITIAL  TEMPERATURE. 

C 

C  SET  THE  INITIAL  CONDITION  OF  THE  REGENERATOR  AND  GAS  AT  A 

C  LINEAR  DIFTRIBUTION  BETWEEN  TOI  AND  TLI. 

TRP (1)=T0I 
TGP (1)=T0I 
TTP (1)=T0I 

C  THE  MATRIX  RUNS  FROM  N=2  TO  NMAX-1 . 

DO  26  N=2, NMAX-1 

TRP (N)  =  (TLI-T0I) /  (NMAX-3) * (N-2)+T0I 
TGP  (N)  =TRP  (N) 

TTP  (N)=TRP  (N) 


no  now  ooooooooo  no  noon  noon  no  noon 


26  CONTINUE 

TRP (NMAX) =TLI 
TGP (NMAX)=TLI 
TTP (NMAX)=TLI 

AT  SOME  EARLY  TIME,  THE  TEMPERATURE  AT  N=1  IS  IMPULSIVELY 
CHANGED  TO  A  WARMER  TEMPERATURE. 

N=1 
1=1 

DETERMINE  THE  REYNOLDS  AND  NUSSELT  NUMBER. 

REP  S=P  S  *DOTM/AF /ENU/RHO 
HNU=HTC*DH/EKG 

WRITE(*,*)  AST, DOTM, REPS, HNU,  '  AST,  DOTM, REPS, HNU' 

WRITE(*,*)  DP,  AF,ENU,RHO,HNU,EKG,  '  DP, AF, ENU, RHO, HNU,  EKG' 

DEFINE  SOME  OTHER  CONSTANTS  NEEDED  BELOW. 

TA=2 . 0*AS*HTC/DOTM/CP 
TB=2 . 0*RHO*AF*C/DOTM/GAM 
TC=HTC*AS/C/CM/RHOM/ (A* (1-POR) ) 

TJ=2 . 0*HTC*AST/DOTM/CP 
TD=TB-TA-TJ 
TE=EKS/DELX/RHOM/CM/C 
TF=HTC*AST/RHOM/CM/ACXT/C 
TH=1 . 0-2 . 0*TE-TF 

WRITE (*, *)  TA,TB,TC,TD,TE,TF, TH,TJ, 'ta,tb,  .  .' 

START  THE  INTERPOLATION  COUNTER,  ICOUNT,  AND  SET  THE  INITIAL 
INTERPOLATION  MARKER,  10. 

I0=II 
IC0UNT=1 

TG2 (I0)=TGP (NMAX) 

DO  31  I=2,ITRAN 

IF (ICOUNT. EQ.IINT)  THEN 
10=10+1 
ICOUNT=0 
END  IF 

INTERPOLATE  THE  VALUE  OF  THE  INCOMING  GAS  BETWEEN  EXPERIMENTALLY 
MEASURED  DATA  POINTS. 

TEMP=  (  (TGO  (lO+l)  -TGO  (10) )  /IINT)  *ICOUNT+TG0  (10) 

SET  THE  INLET  TEMPERATURE  OF  THE  GAS  AT  THE  INLET  TO  THE  MATRIX 
EQUAL  TO  THE  INTERPOLATED  EXPERIMENTAL  TGO. 

TGC(1)=TEMP 

DETERMINE  THE  VALUE  OF  TG,  TR,AND  TT  AT  THE  NEXT  LOCATION. 

DO  32  N=2,NMAX-1 

TGC  (N)  =  (TGP  (N-1)  -TGP  (N+1)  +TA*TRP  (N)  +TJ*TTP  (N)  +TD*TGP  (N)  )  /TB 
TRC  (N)  =TC*TGP  (N)  +TRP  (N)  *  (1 . 0-TC) 

TTC  (N)  =TE*TTP  (N+1)  +TE*TTP  (N-1)  +TF*TGP  (N)  +TH*TTP  (N) 

CONTINUE 

SET  THE  TEMPERATURE  OF  THE  GAS/MATIX  AS  IT  LEAVES. 

TGC (NMAX) =TGC (NMAX-1) 

TRC (NMAX) =TRC (NMAX-1) 

TTC (NMAX) =TTC (NMAX-1) 

PUT  AN  ADIABATIC  CONDITION  AT  THE  FRONT  OF  THE  TUBE  AND  MATRIX. 
TRC(1)=TRC(2) 

TTC(1)=TTC(2) 


oooooo  o  oonocoo  n  on  ooo 


C  SET  THE  CURRENT  TIME  STEP  TEMPERATURES  TO  THE  PREVIOUS  OPES. 

DO  55  J=1,NMAX 
TCP (J)=TGC(J) 

TRP (J)=TRC(J) 

TTP (J) =TTC(J) 

55  CONTINUE 

LET  TG2  EQUAL  THE  CALCULATED  TEMPERATURE  AT  THE  SAME  TIMES 
AS  THE  EXPERIMENTAL  TEMPERATURES  AT  THE  EXIT. 

IF(ICOUNT.EQ.O)  THEN 
TG2 (IO)=TGC(NMAX) 

END  IF 

CHECK  FOR  INSTABILITY. 

IF  (ABS (TG2 (10) ) .GE. 100. )  THEN 

WRITE  (*,*)  '  INSTABILITY  OCCURING' 

STOP 
END  IF 

ICOUNT=ICOUNT+l 
CONTINUE 

FIND  THE  VALUE  OF  THE  SLOPE  AT  I=LOC. 

WRITE  (*,*)  LOG 

CSLOPE= (TG2 (LOC-2) +8 . 0*TG2 (LOC+1) -8 . 0*TG2 (LOC-1) 

$-TG2 (LOC+2) ) /12 . 0/SAMPT 

WRITE (*,*)  CSLOPE, 'CSLOPE*******' 

FUNSMS=CSLOPE-ESLOPE 
WRITE (*,*)  'FUNSMS  =  ',FUNSMS 
RETURN 
END 

★**************************************-********************** 

THIS  SUBROUTINE  FINDS  THE  VALUE  OF  HTC  WHICH  BEST  FITS  THE 
VALUE  OF  THE  EXPERIMENTALLY  MEASURED  SEDT. 

DOUBLE  PRECISION  FUNCTION  SZERDT  (HTCMIN,  HTCMAX,  FUNSDT,  TOL,  12, 
$TAU2) 

DOUBLE  PRECISION  HTCMIN, HTCMAX, FUNSDT, TOL, TAU2 
INTEGER  12 
C 

C  A  ZERO  OF  THE  FUNCTION  F  (X)  IS  COMPUTED  IN  THE  INTERVAL  AX,BX 

C 

C  INPUT 
C 

C  AX 
C  BX 
C  F 
C 

C  TOL 
C 
C 
C 

C  OUTPUT 
C 

C  ZEROIN  ABCISSA  APPROXIMATING  A  ZERO  OF  F  IN  THE  INTERVAL  AX,BX 
C 

C  IT  IS  ASSUMED  THAT  F  (AX)  AND  F  (BX)  HAVE  OPPOSITES  SIGNS 

C  WITHOUT  A  CHECK.  ZEROIN  RETURNS  A  ZERO  X  IN  THE  GIVEN  INTERVAL 

C  AX,BX  TO  WITHIN  A  TOLERANCE  4*MACHEPS*ABS (X)  +  TOL,  WHERE  MACHEPS 
C  IS  THE  RELATIVE  MACHINE  PRECISION. 

C  THIS  FUNCTION  SUBPROGRAM  IS  A  SLIGHTLY  MODIFIED  TRANSLATION  OF 


LEFT  ENDPOINT  OF  INITIAL  INTERVAL 
RIGHT  ENDPOINT  OF  INITIAL  INTERVAL 

FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F  (X)  FOR  ANY  X  IN 
THE  INTERVAL  AX,BX 

DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THE 
FINAL  RESULT 


C  THE  ALGOL  60  PROCEDURE  ZERO  GIVEN  IN  R.  BRENT,  ALGORITHMS  FOR 

C  MINIMIZATION  WITHOUT  DERIVATIVES,  PRENTICE-HALL,  INC.  (1973) 

C 

C 

DOUBLE  PRECISION  A,  B, C, D, E, EPS, FA,FB,FC, TOLl,  XM,  P, Q,  R, S 
C 

EXTERNAL  FUNSDT 
C 

C  COMPUTE  EPS,  THE  RELATIVE  MACHINE  PRECISION 
C 

EPS  =1.0 

61  EPS=EPS/2.0 
TOLl  =  1.0  +  EPS 

IF  (TOLl. GT. 1.0) GO  TO  61 
C 

C  INITIALIZATION 
C 

A=HTCMIN 

B=HTCMAX 

FA=FUNSDT  ( 12 , A, TAU2 ) 

FB=FUNSDT ( 12 , B, TAU2 ) 

WRITE (*,*)  '  fa=  ',FA, '  fb=  ',FB 
C 

IF  (  (FA/FB)  .GT.O .0)  GO  TO  71 
C 

C  BEGIN  STEP 

C 

62  C=A 
FC=FA 
D=B-A 
E=D 

63  IF(ABS(FC) .GE.ABS(FB))GO  TO  64 

A=B 

B=C 

C=A 

FA=FB 

FB=FC 

FC=FA 

C 

C  CONVERGENCE  TEST 
C 

64  TOLl  =  2.0*EPS*ABS(B)  +  0.5*TOL 

XM  =  0.5* (C-B) 

IF(ABS(XM)  .LE.TOLDGO  TO  69 
IF(FB.EQ.0.0)GO  TO  69 
C 

C  IS  BISECTION  NECESSARY? 

C 

IF(ABS(E)  .LT.TOLDGO  TO  67 
IF(ABS(FA) .LE.ABS(FB))GO  TO  67 
C 

C  IS  QUADRATIC  INTERPOLATION  POSSIBLE? 

C 

IF(A.NE.C)GO  TO  65 
C 

C  LINEAR  INTERPOLATION 
C 

S=FB/FA 
P=2.0*XM*S 
Q=1.0-S 
GO  TO  60 
C 

C  INVERSE  QUADRATIC  INTERPOLATION 
C 

65  Q=FA/FC 
R=FB/FC 


S=FB/FA 

P=S* (2 . 0*XM*Q* (Q-R) - (B-A) * (R-1 . 0)  ) 

Q  =  (Q  -  1.0)*(R  -  1.0)*(S  -  1.0) 

C 

C  ADJUST  SIGNS 
C 

60  IF(P.GT.0.0)Q  =  -Q 

P  =  ABS(P) 

C 

C  IS  INTERPOLATION  ACCEPTABLE? 

C 

IF( (2.0*P) .GE. (3.0*XM*Q  -  ABS (T0L1*Q) ) ) GO  TO  67 
IF(P.GE.ABS(0.5*E*Q) )GO  TO  67 
E=D 
D=P/Q 
GO  TO  68 
C 

C  BISECTION 
C 

67  D=XM 
E=D 

C 

C  COMPLETE  STEP 
C 

68  A=B 
FA=FB 

IF  (ABS  (D)  .GT. TOLD  B  =  B  +  D 

IF  (ABS  (D)  .LE. TOLD  B  =  B  +  SIGN (TOLl,  XM) 

FB  =  FUNSDT(I2,B,TAU2) 

WRITE (*,*)  '  The  New  Value  of  fb  is  '  ,B 
IF( (FB*  (FC/ABS(FC) ) ) .GT.O.O)GO  TO  62 
GO  TO  63 
C 

C  DONE 
C 

69  SZERDT  =  B 
RETURN 

71  WRITE (*,*)  'FUNSDT(BAX)/FUNSDT(BBX> .GT.0.0' 

STOP 

END 

C 

C 

DOUBLE  PRECISION  FUNCTION  FUNSDT (12, HTC, TAU2) 

C  THIS  IS  A  FUNCTION  SUBROUTINE  WHICH  DETERMINES  THE  VALUE 

C  OF  THE  DIFFERENCE  BETWEEN  THE  EXPERIMENTAL  AND  CALCULATED 

C  SPONGE  EFFECT  DELAY  TIME,  SEDT,  OF  THE  TIME-TEMPERATURE 

C  CURVE  TRANSIENT. 

C  IT  IS  USED  TO  GET  THE  BEST  VALUE  FOR  HTC  IN  CONJUNCTION 

C  WITH  THE  FUNCTION  SZERDT  (HTCMIN,  HTCMAX,  FA/FB,  TOL)  . 

C  msieSc.f  IS  THE  CONTROLING  PROGRAM  UNIT.  SEE  PROGRAM  tsie.f 

C  FOR  OTHER  DETAILS. 

C 

C  TA, TB,TC, TD,TE,TF,TH,TI  ARE  CONSTANTS  USED  IN  CALCULATIONS. 

C  CSEDT  IS  THE  CALCULATED  SEDT. 

C  SEDT  IS  THE  EXPERIMENTALLY  MEASURED  SEDT. 

C  TAU2  IS  THE  TIME  AT  THE  END  OF  THE  SEDT. 

C  12  IS  THE  INDEX  FOR  TAU2;  I2C  IS  THE  CALCULATED  12. 

C 

C 

DOUBLE  PRECISION  TGC (155) , TRC (155) , TTC (155) ,  TGP  (155) , TRP (155) , 
$TTP  (155)  ,  TGO (20001) , TGL  (20001) , TG2 (20001) 

INTEGER  I,N,NMAX, II, IF, 12, 10, ITRAN, IINT, ICOUNT,  LOC, I2C,  J,  K 
DOUBLE  PRECISION  RHO, TA, TB, TC, TD, CM,  RHOM, A,  CP, GAM, 

$DELT,  DELX, DH, AF, ENU, POR, D, PI, DP, PS, EKG, EKS, 

$TE, TF,  TH, TJ, TOI, TLI, SAMPT,  TG2AV, 


ooono  oo  onon  noon  oo  oono 


$AS,  C,  HTC,  DOTH,  REPS,  TEMP,  THK,  ACXT,  AST,  SEDT,  CSEDT,  TAU2 
C 

COMMON  TGO,  TGL,  DOTM,  AS,  CP,  DELX,  SAMPT,  RHO,  AT,  C,  CM,  GAM,  DELT,  THK, 
$POR,RHOM,  TOI,TLI,DP,ENU,EKG,DH,PS,  II,  IF,  IINT,  ITRAN, 

$D, PI, A, ACXT, AST,  SKS, NMAX, LOC 
C 

INTRINSIC  COS,MAX, SQRT,ABS,NINT,ATAN 
C  AT  THE  INITIAL  TIME,  THE  REGENERATOR  AND  GAS  ARE  AT  A 

C  STEADY  STATE  INITIAL  TEMPERATURE. 

C 

C  SET  THE  INITIAL  CONDITION  OF  THE  REGENERATOR  AND  GAS  AT  A 

C  LINEAR  DIFTRIBUTION  BETWEEN  TOI  AND  TLI . 

TRP (1)=T0I 
TGP (1)=T0I 
TTP (1)=T0I 

C  THE  MATRIX  RUNS  FROM  N=2  TO  NMAX-1. 

DO  26  N=2, NMAX-1 

TRP (N)= (TLI-TOI) / (NMAX-3) * (N-2)+T0I 
TGP (N)=TRP (N) 

TTP  (N)  =TRP  (N) 

2  6  CONTINUE 

TRP (NMAX) =TLI 
TGP (NMAX) =TLI 
TTP (NMAX)=TLI 

AT  SOME  EARLY  TIME,  THE  TEMPERATURE  AT  N=1  IS  IMPULSIVELY 
CHANGED  TO  A  WARMER  TEMPERATURE. 


N=1 

1=1 


DETERMINE  THE  REYNOLDS  AND  NUSSELT  NUMBER. 

REPS=PS*DOTM/AF/ENU/RHO 

HNU=HTC*DH/EKG 

WRITE(*,*)  AST,  DOTM,  REPS,  HNU,  '  AST,  DOTM,  REPS,  HNU' 

WRITE  (*,*)  DP,AF,ENU,RHO,HNU,EKG,  '  DP,  AF,  ENU,  RHO,  HNU,  EKC' 

DEFINE  SOME  OTHER  CONSTANTS  NEEDED  BELOW. 

TA=2 . 0*AS*HTC/DOTM/CP 
TB=2 . 0*RHO*AF*C/DOTM/GAM 
TC=HTC*AS/C/CM/RHOM/ (A* (1-POR)  ) 

TJ=2 . 0*HTC*AST/DOTM/CP 

TD=TB-TA-TJ 

TE=EKS/DELX/RHOM/CM/C 

TF=HTC*AST/RHOM/CM/ACXT/C 

TH=1.0-2.0*TE-TF 

WRITE  (*,  *)  TA,  TB,  TC,TD,TE,TF,TH,TJ, 'ta,  tb,  .  .' 

START  THE  INTERPOLATION  COUNTER,  ICOUNT,  AND  SET  THE  INITIAL 
INTERPOLATION  MARKER,  10. 

I0=II 

ICOUNT=l 

TG2 (I0)=TGP (NMAX) 

DO  31  1=2, ITRAN 

IF (ICOUNT. EQ. IINT)  THEN 
10=10+1 
ICOUNT=0 
END  IF 

INTERPOLATE  THE  VALUE  OF  THE  INCOMING  GAS  BETWEEN  EXPERIMENTALLY 
MEASURED  DATA  POINTS. 

TEMP=(  (TGO  (I0+1)-TG0  (10)  )  /IINT)  *ICOUNT+TGO  (10) 

SET  THE  INLET  TEMPERATURE  OF  THE  GAS  AT  THE  INLET  TO  THE  MATRIX 


C  EQUAL  TO  THE  INTERPOLATED  EXPERIMENTAL  TGO . 

TGC(1)=TEMP 

C 

C  DETERMINE  THE  VALUE  OF  TG,  TR,AND  TT  AT  THE  NEXT  LOCATION. 

DO  32  N=2,NMAX-1 
C 

TGC  (N)  =  (TD*TGP  (N)  +TGP  (N-l)-TGP  (N+1)  +TA*TRP  (N)  +TJ*TTP  (N)  )  /TB 
TRC  (N)  =TC*TGP  (N)  +TRP  (N)  *  (1 . 0-TC) 

TTC  (N)  =TE*TTP  (N+1)  +TE*TTP  (N-1)  +TF*TGP  (N)  +TH*TTP  (N) 

32  CONTINUE 

C 

C  SET  THE  TEMPERATURE  OF  THE  GAS/MATIX  AS  IT  LEAVES. 

TGC (NMAX) =TGC (NMAX-1) 

TRC (NMAX) =TGC (NMAX-1) 

TTC (NMAX) =TTC (NMAX-1) 

C 

C  PUT  AN  ADIABATIC  CONDITION  AT  THE  FRONT  OF  THE  TUBE  AND  MATRIX. 

TRC(1)=TRC(2) 

TTC(1)=TTC(2) 

C 

C  SET  THE  CURRENT  TIME  STEP  TEMPERATURES  TO  THE  PREVIOUS  OPES. 

DO  55  J=1,NMAX 
TGP (J)=TGC(J) 

TRP (J)=TRC(J) 

TTP (J)=TTC(J) 

55  CONTINUE 

C 

C  LET  TG2  EQUAL  THE  CALCULATED  TEMPERATURE  AT  THE  SAME  TIMES 

C  AS  THE  EXPERIMENTAL  TEMPERATURES  AT  THE  EXIT. 

IF (ICOUNT.EQ.O)  THEN 
TG2  (I0)=TGC (NMAX) 

END  IF 
C 

C  CHECK  FOR  INSTABILITY. 

IF(ABS(TG2(I0)) .GE.IOO.)  THEN 

WRITE (*,*)  '  INSTABILITY  OCCURRING' 

STOP 
END  IF 
C 

ICOUNT=ICOUNT+l 

C 

31  CONTINUE 

C 

C  FIND  THE  VALUE  OF  THE  CSEDT. 

DO  85  J=II,IF 
TG2AV=0 . 0 

DO  87  K=II,J 
TG2AV=TG2AV+TG2 (K) 

87  CONTINUE 

TG2AV=TG2AV/ ( J-II+1) 

IF  ( (TG2 (J+50)-TG2AV) .GT.0.4)  THEN 

CSEDT= (J-II) *SAMPT 

I2C=J 

GO  TO  86 

END  IF 

85  CONTINUE 

86  CONTINUE 
C 

C  DETERMINE  THE  VALUE  OF  SEDT. 

SEDT= (I2-II) *SAMPT 
C 

FUNSDT=I2C-I2 

WRITE (*,*)  'FUNSDT  =  ',FUNSDT 

WRITE (*,*)  'CSEDT  =  ', CSEDT, '  SEDT=  ' , SEDT 

RETURN 

END 


cr>000  O  noo  m  OOO 


C 

c 
c 

C  THIS  SUBROUTINE  FINDS  THE  VALUE  OF  HTC  WHICH  MINIMIZES  THE 

C  ROOT  MEAN  SQUARE  DIFFERENCE  BETWEEN  THE  CALCULATED  AND 

C  EXPERIMENTAL  VALUES  OF  THE  OUTLET  TEMPERATURE  AFTER  SEDT. 

C 

DOUBLE  PRECISION  FUNCTION  SZERRM (HTCMIN, HTCMAX,  FUNSRM, TOL,  12 ) 
DOUBLE  PRECISION  HTCMIN, HTCMAX,  FUNSRM,  TOL 
INTEGER  12 
C 

C  A  ZERO  OF  THE  FUNCTION  F  (X)  IS  COMPUTED  IN  THE  INTERVAL  AX,BX 
C 

C  INPUT 
C 

C  AX  LEFT  ENDPOINT  OF  INITIAL  INTERVAL 
C  BX  RIGHT  ENDPOINT  OF  INITIAL  INTERVAL 

C  F  FUNCTION  SUBPROGRAM  WHICH  EVALUATES  F  (X)  FOR  ANY  X  IN 

C  THE  INTERVAL  AX, BX 

C  TOL  DESIRED  LENGTH  OF  THE  INTERVAL  OF  UNCERTAINTY  OF  THE 
C  FINAL  RESULT 

C 
C 

C  OUTPUT 
C 

C  ZEROIN  ABCISSA  APPROXIMATING  A  ZERO  OF  F  IN  THE  INTERVAL  AX,BX 
C 

C  IT  IS  ASSUMED  THAT  F  (AX)  AND  F  (BX)  HAVE  OPPOSITES  SIGNS 

C  WITHOUT  A  CHECK.  ZEROIN  RETURNS  A  ZERO  X  IN  THE  GIVEN  INTERVAL 

C  AX,BX  TO  WITHIN  A  TOLERANCE  4*MACHEPS*ABS (X)  +  TOL,  WHERE  MACHEPS 
C  IS  THE  RELATIVE  MACHINE  PRECISION. 

C  THIS  FUNCTION  SUBPROGRAM  IS  A  SLIGHTLY  MODIFIED  TRANSLATION  OF 
C  THE  ALGOL  60  PROCEDURE  ZERO  GIVEN  IN  R.  BRENT,  ALGORITHMS  FOR 

C  MINIMIZATION  WITHOUT  DERIVATIVES,  PRENTICE-HALL,  INC.  (1973) 

C 

C 

DOUBLE  PRECISION  A, B, C, D, E,EPS, FA, FB, FC, TOLL, XM, P, Q, R, S 
C 

EXTERNAL  FUNSRM 

COMPUTE  EPS,  THE  RELATIVE  MACHINE  PRECISION 

EPS  =1.0 
EPS=EPS/2.0 
TOLl  =  1.0  +  EPS 
IF (TOLl.GT.l.O)GO  TO  61 

INITIALIZATION 

A=HTCMIN 
B=HTCMAX 
FA=FUNSRM(I2,A) 

FB=FUNSRM(I2,B) 

WRITE (*,*)  '  fa=  ',FA, '  fb=  ',FB 

IF( (FA/FB) .GT.0.0)  GO  TO  71 

BEGIN  STEP 

C=A 
FC=FA 
D=B-A 
E=D 

IF(ABS(FC) .GE.ABS(FB))GO  TO  64 
A=B 
B=C 


63 


C=A 

FA=FB 

FB=FC 

FC=FA 

C 

C  CONVERGENCE  TEST 
C 

64  TOLl  =  2.0*EPS*ABS(B)  +  0.5*TOL 
XM  =  0.5* (C-B) 

IF(ABS  (XM)  .LE.TOLDGO  TO  69 
IF(FB.EQ.O.O)GO  TO  69 
C 

C  IS  BISECTION  NECESSARY? 

C 

IF(ABS(E)  .LT.TOLDGO  TO  67 
IF (ABS (FA) .LE.ABS (FB) )GO  TO  67 
C 

C  IS  QUADRATIC  INTERPOLATION  POSSIBLE? 

C 

IF(A.NE.C)GO  TO  65 
C 

C  LINEAR  INTERPOLATION 
C 

S=FB/FA 
P=2.0*XM*S 
Q=1.0-S 
GO  TO  60 
C 

C  INVERSE  QUADRATIC  INTERPOLATION 
C 

65  Q=FA/FC 
R=FB/FC 
S=FB/FA 

P=S* (2 . 0*XM*Q* (Q-R) - (B-A) * (R-1 . 0) ) 

Q  =  (Q  -  1.0)*(R  -  1.0)*(S  -  1.0) 

C 

C  ADJUST  SIGNS 
C 

60  IF (P.GT.O.O)Q  =  -Q 

P  =  ABS(P) 

C 

C  IS  INTERPOLATION  ACCEPTABLE? 

C 

IF( (2.0*P) .GE. (3.0*XM*Q  -  ABS (TOLl*Q) ) ) GO  TO  67 
IF(P.GE.ABS(0.5*E*Q) )GO  TO  67 
E=D 
D=P/Q 
GO  TO  68 
C 

C  BISECTION 
C 

67  D=XM 
E=D 

C 

C  COMPLETE  STEP 
C 

68  A=B 
FA=FB 

IF (ABS (D) .GT. TOLl)  B  =  B  +  D 

IF (ABS (D)  .LE. TOLl)  B  =  B  +  SIGN (TOLl, XM) 

FB  =  FUNSRM(I2,B) 

WRITE (*,*)  '  The  New  Value  of  B  is  ',B 
IF( (FB*  (FC/ABS(FC) ) )  .GT.0.0)GO  TO  62 
GO  TO  63 
C 

C  DONE 


c 

69  SZEKRM  =  B 

RETURN 

71  WRITE (*,*)  'FUNSRM(BAX) /FUNSRM(BBX) .GT.0.0' 

STOP 

END 

C 

C 

DOUBLE  PRECISION  FUNCTION  FUNSRM (12,  HTC) 

C  THIS  IS  A  FUNCTION  SUBROUTINE  WHICH  DETERMINES  THE  VALUE 

C  OF  THE  ROOT  MEAN  SQUARE  DIFFERENCE  BETWEEN  THE  EXPERIMENTAL 

C  AND  CALCULATED  VALUES  OF  THE  OUTLET  TEMPERATURE 

C  OF  THE  TIME-TEMPERATURE  CURVE  TRANSIENT. 

C  IT  IS  USED  TO  GET  THE  BEST  VALUE  FOR  HTC  IN  CONJUNCTION 

C  WITH  THE  FUNCTION  SZERRM(HTCMIN,  HTCMAX,  FA/FB,  TOL)  . 

C  msieSc.f  IS  THE  CONTROLING  PROGRAM  UNIT.  SEE  PROGRAM  tsie.f 

C  FOR  OTHER  DETAILS. 

C 

C  TA,  TB,  TC,  TD,  TE,  TF,  TH,  TI  ARE  CONSTANTS  USED  IN  CALCULATIONS. 

C  12  IS  THE  INDEX  OF  THE  TIME  AT  THE  END  OF  THE  SEDT . 

C  RMSD/1/2  ABE  ROOT  MEAN  SQUARED  DIFFERENCES. 

C 

C 

DOUBLE  PRECISION  TGC (155) , TRC (155) , TTC (155)  ,  TGP (155) , TRP  (155) , 
$TTP  (155) , TGO (20001) , TGL (20001) , TG2  (20001) 

INTEGER  I,N,NMAX,  II,  IF,  12, 10,  ITRAN,  IINT,  ICOUNT,  LOHI,  LOC 
DOUBLE  PRECISION  RHO, TA, TB, TC, TD, CM, RHOM,  A,  CP,  GAM, 

$DELT,  DELX,  DH,  AF,  ENU,  POR,  D,  PI,  DP,  PS,  EKG,  EKS, 

$TE, TF,  TH, TJ, TOI, TLI, SAMPT, 

$AS,  C,  HTC,DOTM,  REPS,  TEMP,  THK,  ACXT,  AST,RMSD,  RMSDl,  RMSD2,  DT 
C 

COMMON  TGO,  TGL,  DOTM,  AS,  CP,  DELX,  SAMPT,  RHO,  AF,  C,  CM,  GAM,  DELT,  THK, 
$POR,  RHOM,  TOI,  TLI, DP ,  ENU, EKG, DH, PS,  II,  IF,  IINT,  ITRAN, 

$D , P I ,  A ,  ACXT , AST ,  EKS , NMAX,  LOC 
C 

INTRINSIC  COS,MAX, SQRT,ABS,NINT,ATAN 
C 

C  FIND  THE  SLOPES  OF  THE  RMSD  VERSUS  HTC  CURVE  AT  HTC  AND 

C  HTC  +  10.0. 

C 

LOHI=l 
82  N=1 

1=1 

C  AT  THE  INITIAL  TIME,  THE  REGENERATOR  AND  GAS  ARE  AT  A 

C  STEADY  STATE  INITIAL  TEMPERATURE. 

C 

C  SET  THE  INITIAL  CONDITION  OF  THE  REGENERATOR  AND  GAS  AT  A 

C  LINEAR  DIFTRIBUTION  BETWEEN  TOI  AND  TLI. 

TRP (1)=T0I 
TGP (1)=T0I 
TTP (1)=T0I 

C  THE  MATRIX  RUNS  FROM  N=2  TO  NMAX-1. 

DO  26  N=2, NMAX-1 

TRP  (N)  =  (TLI-TOI)  /  (NMAX-3)  *  (N-2)  +T0I 
TGP  (N)  =TRP  (N) 

TTP  (N)  =TRP  (N) 

26  CONTINUE 

TRP  (NMAX)  =TLI 
TGP (NMAX) =TLI 
TTP (NMAX) =TLI 
C 

C  AT  SOME  EARLY  TIME,  THE  TEMPERATURE  AT  N=1  IS  IMPULSIVELY 

C  CHANGED  TO  A  WARMER  TEMPERATURE. 

C 

C  DETERMINE  THE  REYNOLDS  AND  NUSSELT  NUMBER. 


REP  S=P  S  *DOTM/AE /ENU/RHO 
HNU=HTC*DH/EKG 

C  WRITE  (*,*)  AST, DOTH, REPS, HNU,  '  AST,  DOTH, REPS,  HNU' 

C  WRITE(*,*)  DP,AF,ENU,RHO,HNU,EKG, '  DP, AF, ENU, RHO, HNU, EKG' 

C 

C  DEFINE  SOME  OTHER  CONSTANTS  NEEDED  BELOW. 

TA=2 . 0*AS*HTC/DOTM/CP 
TB=2 . 0*RHO*AF*C/DOTM/GAM 
TC=HTC*AS/C/CM/RHOM/ (A* (1-POR) ) 

TJ=2 . 0*HTC*AST/DOTM/CP 

TD=TB-TA-TJ 

TE=EKS/DELX/RHOM/CM/C 

TF=HTC*AST/RHOM/CM/ACXT/C 

TH=1.0-2.0*TE-TF 

C  WRITE (*,*)  TA, TB,TC,TD,TE,TF,  TH,TJ, 'ta,tb,  . 

C 

C  START  THE  INTERPOLATION  COUNTER,  ICOUNT,  AND  SET  THE  INITIAL 

C  INTERPOLATION  MARKER,  10. 

I0=II 

ICOUNT=l 

TG2 (I0)=TGP (NMAX) 

C 

DO  31  1=2, ITRAN 
C 

IF(ICOUNT.EQ.IINT)  THEN 
10=10+1 
ICOUNT=0 
END  IF 
C 

C  INTERPOLATE  THE  VALUE  OF  THE  INCOMING  GAS  BETWEEN  EXPERIMENTALLY 

C  MEASURED  DATA  POINTS. 

TEMP= ( (TGO (lO+l) -TGO (10) ) /IINT) *ICOUNT+TG0 (10) 

C 

C  SET  THE  INLET  TEMPERATURE  OF  THE  GAS  AT  THE  INLET  TO  THE  MATRIX 

C  EQUAL  TO  THE  INTERPOLATED  EXPERIMENTAL  TGO. 

TGC(1)=TEMP 

C 

C  DETERMINE  THE  VALUE  OF  TG,TR,AND  TT  AT  THE  NEXT  LOCATION. 

DO  32  N=2,NMAX-1 
C 

TGC  (N)  =  (TD*TGP  (N)  +TGP  (N-1)  -TGP  (N+1)  +TA*TRP  (N)  +TJ*TTP  (N)  )  /TB 
TRC (N) =TC*TGP (N) +TRP (N) * (1 . 0-TC) 

TTC  (N)  =TE*TTP  (N+1)  +TE*TTP  (N-1)  +TF*TGP  (N)  +TH*TTP  (N) 

32  CONTINUE 

C 

C  SET  THE  TEMPERATURE  OF  THE  GAS/MATIX  AS  IT  LEAVES. 

TGC (NMAX) =TGC (NMAX-1 ) 

TRC  (NMAX)  =TGC  (NMAX-1 ) 

TTC (NMAX) =TTC (NMAX-1 ) 

C 

C  PUT  AN  ADIABATIC  CONDITION  AT  THE  FRONT  OF  THE  TUBE  AND  MATRIX. 

TRC(1)=TRC(2) 

TTC(1)=TTC(2) 

C 

C  SET  THE  CURRENT  TIME  STEP  TEMPERATURES  TO  THE  PREVIOUS  OPES. 

DO  55  J=1,NMAX 
TGP (J)=TGC(J) 

TRP ( J) =TRC  ( J) 

TTP (J)=TTC(J) 

55  CONTINUE 

C 

C  LET  TG2  EQUAL  THE  CALCULATED  TEMPERATURE  AT  THE  SAME  TIMES 

C  AS  THE  EXPERIMENTAL  TEMPERATURES  AT  THE  EXIT. 

IF ( ICOUNT. EQ.O)  THEN 
TG2 (I0)=TGC (NMAX) 

END  IF 


OOO  OOoo  oocoon  oo 


CHECK  FOR  INSTABILITY. 

IF (ABS (TG2 (10) ) .GE. 100.)  THEN 

WRITE (*,*)  '  INSTABILITY  OCCURRING' 

STOP 
END  IF 

ICOUNT=ICOUNT+1 
CONTINUE 

DETERMINE  THE  ROOT  MEAN  SQUARED  DIFFERENCE. 

DT  =  0.0 

DO  81  J=I2, (LOC+50) 

DT=(TGL(J)-TG2 (J) ) **2+DT 
CONTINUE 

RMSD=SQRT (DT/ (IF-I2+1) ) 

RECORD  THE  VALUE  OF  RMSD  AND  RUN  FOR  ANOTHER  HTC. 

IF(LOHI.EQ.l)  THEN 
RMSD1=RMSD 
LOHI=2 
HTC=HTC+10. 

GO  TO  82 
ELSE 

RMSD2=RMSD 
END  IF 

DEFINE  FUNSRM  AS  THE  LOCAL  SIMPLE  DIFFERENCE  DERIVATIVE. 

FUNSRM  =  (RMSD2-RMSD1)/10. 

C 

WRITE(*,*)  'FUNSRM  =  ', FUNSRM,'  RMSD1,RMSD2  ARE  '  ,RMSD1,RMSD2 

RETURN 

END 


Appendix  D:  Calibrations 


This  section  describes  the  calibration  of  the  three  most 
important  experimental  measurements :  1 )  temperatures ,  2 ) 
pressure  drops,  and  3)  mass  flow  rate. 

Temperature  Calibration 

The  thermistors  at  the  inlet  and  outlet  were  calibrated 
in-situ  with  a  small  iron-constantan  thermocouple.  The 
thermocouple  was  mounted  midway  in  a  three  inch  empty  tube 
which  was  being  supported  at  each  end  with  the  instrumentation 
connectors.  The  cooler  or  heater  was  turned  on,  and  the 
system  was  allowed  to  come  to  equilibrium  at  seven  different 
temperatures  in  the  range  5-35°C.  A  Labview  calibration  VI 
(My_Thermistor_Calibration. VI)  was  designed  to  read  the  value 
of  the  thermocouple  temperature  and  the  voltage  across  the 
thermistor  bridge,  and  save  this  data  to  a  file.  Some 
hysteresis  was  noted  in  going  from  cold-to-hot  and  then  hot- 
to-cold,  so  a  one  directional  calibration  was  performed.  It 
was  done  going  from  cold  temperatures  to  hot  temperatures, 
since  the  actual  tests  all  had  a  transient  from  cold-to-hot. 

The  collection  of  voltages  and  related  temperatures  were 
saved  to  a  data  file  and  a  least-square  curve  fit  was  done. 
The  calibration  data  were  inserted  into  an  equation  (Doebelin, 
1983:  609-613)  which  has  been  shown  to  describe  accurately  the 
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behavior  of  thermistors  over  temperature  ranges  less  than 
50°C: 


1  ■  0 
T 


A^ilnRj.)^  +  A^(lnR^)^  +  A^ilnR^)  +  A^ 


(30) 


where  is  the  resistance  of  the  thermistor  (Ohms) 

A^  are  constants  determined  by  the 

calibration 
Ai  =  2.59  X  10“^ 

A2  =  -1.94  X  10'^ 

A3  =  7.74  X  10-^ 

A4  =  2.24  X  10'^ 

and  T  is  the  temperature  (°C) 

This  equation  was  used  in  conjunction  with  the  equation 
describing  the  voltage  drop  for  one  leg  of  a  Wheatstone 
bridge.  A  curve  fit  was  made  over  the  range  5°C  to  35°C  for 
each  of  the  two  thermistors.  The  resulting  curve  matched  the 
data  to  ±  0.2°C. 

Pressure  Drop 

The  three  Validyne  differential  pressure  transducers  were 
calibrated  with  a  dead  weight  tester.  A  known  pressure  was 
imposed  across  the  diaphragm  of  each  transducer  and  the 
resulting  voltage  output  was  recorded.  The  calibration  showed 
a  linear  relation  between  voltage  output  and  pressure.  The 
voltage  range  for  each  transducer  was  0-10  VDC.  A  least- 
squares  curve  fit  was  performed  for  each  transducer  and  the 
coefficients  were  incorporated  into  the  Labview  VI  which 
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monitors  the  tests.  These  linear  curve  fit  equations  matched 
the  transducer  calibration  data  to  less  a  standard  deviation 
of  than  1%  of  full  range  for  all  three  cases. 

Mass  Flow  Rate 

The  mass  flow  rate  was  determined  by  a  relationship  which 
includes  the  measured  pressure  drop  across  the  venturi  tube 
constriction,  the  density  of  the  gas,  and  some 
parameters  of  the  geometry  of  the  venturi  tube  (Benedict, 
1984:  Chapter  21)  .  The  accuracy  of  the  relation  for  mass  flow 
rate  was  checked  with  an  independent  experiment.  A  stagnation 
chamber  with  a  nozzle  outlet  was  connected  in-line  with  the 
venturi  tube.  The  gas  flow  in  the  stagnation  chamber  was 
idealized  as  having  zero  velocity,  and  gasdynamic  relations 
for  nozzle  flow  from  stagnation  to  atmospheric  were  used  to 
determine  the  mass  flow  rate  out  of  the  cell  at  equilibrium. 
A  Labview  VI  was  devised  to  compare  the  venturi  tube  results 
to  those  from  the  pressure  cell.  If  a  nozzle  efficiency  of 
0.95  is  used  for  the  stagnation  chamber  outlet,  the  venturi 
tube  arrangement  measured  the  flow  rate  to  within  ±  7.0%  at 
the  lowest  flow  rates,  and  to  within  ±  4.0%  at  the  highest 
Reynolds  number . 
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Appendix  E:  Test  Data 


In  this  appendix,  the  final  values  of  the  Reynolds 
number,  Nusselt  number,  friction  factor,  compactness  factor, 
reduction  factor,  and  effectiveness  for  the  twenty  eight  test 
runs  are  included. 
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Re 

Nu 

f 

CF 

RF 

SPEFF 

7.6543 

0.8010 

1.6281 

0.0830 

1.0000 

0.9460 

14.6904 

1.1218 

1.1638 

0.0847 

1.0000 

0.8820 

8.2452 

0.6672 

1.4964 

0.0698 

0.8465 

0.8730 

15.8082 

0.7285 

1.1426 

0.0521 

0.8465 

0.8820 

8.9725 

0.9508 

1.6349 

0.0837 

0.7456 

0.8990 

9.1630 

1.0130 

1.6179 

0.0882 

0.7456 

0.9060 

16.6818 

1.3136 

1.1780 

0.0863 

0.7456 

0.9480 

15.1212 

1.0475 

2.5664 

0.0348 

0.5526 

0.8190 

27.8053 

0.9822 

2.0428 

0.0223 

0.5526 

0.6980 

39.1208 

2.4477 

0.7860 

0.1027 

1.0000 

0.8340 

75.0169 

3.3876 

0.6623 

0.0880 

1.0000 

0.8550 

40.4648 

2.2300 

0.8273 

0.0860 

0.8574 

0.7940 

74.0514 

2.5543 

0.7299 

0.0610 

0.8574 

0.7730 

45.3124 

2.0922 

0.8700 

0.0685 

0.7013 

0.7490 

47.8531 

2.3574 

0.8449 

0.0753 

0.7013 

0.7640 

91.5990 

3.0476 

0.7507 

0.0572 

0.7013 

0.7790 

75.5909 

3.2338 

2.0015 

0.0276 

0.5000 

0.7150 

86.2827 

3.4506 

1.9613 

0.0263 

0.5000 

0.7120 

15.6136 

1.1117 

0.8838 

0.1040 

1.0000 

0.8790 

15.0196 

1.2263 

0.9003 

0.1170 

1.0000 

0.8980 

37.8461 

2.8948 

0.6658 

0.1483 

1.0000- 

0.9010 

16.8837 

1.0339 

0.9857 

0.0802 

0.8302 

0.8450 

40.9000 

1.7855 

0.7321 

0.0770 

0.8302 

0.9280 

17.8880 

1.2749 

1.1168 

0.0824 

0.7207 

0.8620 

43.7263 

1.7756 

0.8548 

0.0613 

0.7207 

0.7770 

30.4596 

2.0455 

1.8725 

0.0463 

0.5328 

0.8280 

30.8783 

2.0618 

1.8439 

0.0467 

0.5328 

0.8370 

48.8819 

2.9644 

1.6604 

0.0471 

0.5328 

0.8860 
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